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ABSTRACT 

We carry out simulations of gravitationally unstable discs using a Smoothed Particle Hy- 
drodynamics (SPH) code and a grid-based hydrodynamics code, FARGO, to understand the 
previous non-convergent results reported by Meru & Bate ( 201 la| l. We obtain evidence that 



convergence with increasing resolution occurs with both SPH and FARGO and in both cases 
we find that the critical cooling timescale is larger than previously thought. We show that SPH 
has a first-order convergence rate while FARGO converges with a second-order rate. We show 
that the convergence of the critical cooling timescale for fragmentation depends largely on the 
numerical viscosity employed in both SPH and FARGO. With SPH, particle velocity disper- 
sion may also play a role. We show that reducing the dissipation from the numerical viscosity 
leads to larger values of the critical cooling time at a given resolution. For SPH, we find that 
the effect of the dissipation due to the numerical viscosity is somewhat larger than had previ- 
ously been appreciated. In particular, we show that using a quadratic term in the SPH artificial 
viscosity (/?sph) that is too low appears to lead to excess dissipation in gravitationally unsta- 
ble discs, which may affect any results that sensitively depend on the thermodynamics, such 
as disc fragmentation. We show that the two codes converge to values of the critical cooling 
timescale, /3 cr jt > 20 (for a ratio of specific heats of 7 = 5/3), and perhaps even as large as 
/3 crit sa 30. These are approximately 3 — 5 times larger than has been found by most previous 
studies. This is equivalent to a maximum gravitational stress that a disc can withstand without 
fragmenting of aGi.crit ~ 0.013 — 0.02, which is much smaller than the values typically 
used in the literature. It is therefore easier for self-gravitating discs to fragment than has been 
concluded from most past studies. 

Key words: accretion, accretion discs - protoplanetary discs - planets and satellites: formation 
- gravitation - instabilities - hydrodynamics 



1 INTRODUCTION 

Historically, there have been two key quantities that have been used 
to determine whether a self-gravitating disc is likely to fragment. 
The first is the stability parameter (Toomre 19641, 



Gammie (2001 1 showed that in addition to the stability crite- 
rion above, the disc must cool at a fast enough rate. Using shearing 
sheet simulations, he showed that if the cooling timescale can be 
parametrized as 



Q = = (i) 

7TZjG 

where c s is the sound speed in the disc, fi: cp is the epicyclic fre- 
quency, which for Keplerian discs is approximately equal to the 
angular frequency, £7, £ is the surface mass density and G is the 
gravitational constant. Toomre ( 1964) showed that for an infinites- 
imally thin disc to fragment, the stability parameter must be less 
than a critical value, Q cr it ~ 1. 
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where 



tcool — u 



du c 



dt 



(2) 



(3) 



u is the specific internal energy and du coo i/dt is the total spe- 
cific cooling rate, then for fragmentation we require (3 < /3 C rit- 
According to |Gammie| {2001|>, /3 cr j t ~ 3 for a ratio of specific 
heats, 7 — 2. Rice, Lodato & Armitage (2005) carried out three- 
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dimensional simulations using a Smoothed Particle Hydrodynam- 
ics (SPH) code and showed that this cooling parameter is dependent 
on the equation of state. They showed that /3 cr it ~ 6 — 7 for discs 
with 7 = 5/3 and f3 CI a « 12 — 13 for discs with 7 = 7/5. 

|Gammie| ( |2001| > and |Rice et al.| (2005l also showed that in a 
steady state disc where the dominant form of heating is that due to 
gravitational instabilities, since the gravitational stress in a disc can 
be linked to the cooling timescale in the disc using 

«GI = r —, 7T "J , (4) 

9 7(7- 1) 

the rapid cooling required for fragmentation, p CT n, can also be in- 
terpreted as a maximum gravitational stress that a disc can support 
without fragmenting, which they showed to be aGI,crit ~ 0.06. 

Recently, Meru & Bate (2011a) showed using SPH calcula- 
tions of gravitationally unstable discs similar to those that have 
been performed by |Rice et al.| ([20051 that the previous results on 
the critical cooling timescale had not converged. In particular, they 
found that the critical value of the cooling timescale, P CT n , below 
which a disc would fragment increased linearly with increasing 
spatial resolution. This implied that the critical cooling rate might 
be much greater than that found from past studies (which would, 
for example, have implications for where in a real disc planets may 
form by the gravitational instability method). It also opened the 
question of whether or not a critical cooling rate indeed exists. In- 
stead, a self-gravitating disc that is subject to a fixed cooling rate 
might fragment regardless of the value, given sufficient resolution 
(i.e. a disc may never be able to settle into a self-regulated state). 

|Lodato & Clarke] |20TT| suggested that the non-convergent 
results may be an artefact of SPH artificially smoothing the den- 
sity enhancements, or may be due to artificial viscosity if its ef- 
fect was much larger than expected. |Rice e t al. (2012) suggested 
that the implementation of cooling in SPH may be to blame for 
the lack of convergence. However, Paardekooper, Baruteau & Meru 
pOTT) showed using the two-dimensional grid-based hydrodynam- 
ics code, FARGO, that the non-convergent problem was not specific 
to SPH. The source of non-convergence therefore cannot be con- 
strained to SPH or to three-dimensional codes. Paardekoope r et akj 
< |20 1 1 ) > suggested that the boundary between the turbulent inner disc 
region and the smooth outer disc region (a consequence of starting 
the simulations with smooth initial conditions) may be the cause of 
the non-convergent results presented by Meru & Bate (201 la| >. 

|Bate| p01 1) carried out SPH simulations of the collapse of 
molecular clouds to form protostars and discs. For particular ini- 
tial conditions that lead to disc fragmentation, he noted that higher 
resolution simulations resulted in more fragments. Unlike the simu- 
lations of gravitationally unstable isolated protoplanetary discs dis- 
cussed above, these simulations were of very early stage discs that 
formed prior to stellar core formation and were subject to rapid 
accretion from the surrounding molecular envelope. However, the 
interesting aspect here is that the fragmentation is more prevalent 
in higher resolution simulations of discs modelled using both iso- 
lated discs as initial conditions and using a parametrized cooling 
function (Meru & Bate 2011a) as well as discs formed in molec- 
ular cloud collapse simulations using radiative transfer where no 
such smooth initial conditions are involved (Bate [20lT) . 

|Meru & Bate (201 la) expressed a concern about a lack of con- 
vergence with numerical resolution. However, even if convergence 
with numerical resolution is achieved, convergence between differ- 
ent numerical models is also important, i.e. the result is not believ- 
able if two different codes that can, in principle, model the same 
physical processes, produce physically different results. 



In this paper, we present additional SPH results to those pre- 
sented by Meru & Bate (2011a). Rather than confining our inves- 
tigations to a single hydrodynamics code, we also carry out a code 
comparison by performing further calculations using the grid-based 
Eulerian hydrodynamics code, FARGO. We particularly focus on 
the dependence of the critical cooling timescale on the artificial 
viscosities employed in both codes. 

In Section [2] we describe the numerical methods adopted and 
discuss how numerical viscosity may affect the critical cooling 
timescale in discs in Section [3] We describe the simulations per- 
formed and present our results in Sections H] and B] respectively. 
We discuss and make conclusions in Sections|6]and|7| respectively. 



2 NUMERICAL METHODS 

Our SPH simulations are carried out using the exact same code as 
that used by |Meru & Bate| (2011a|>, originally dev eloped b y |Benz| 
( 1990 ), further developed by |Bate, Bonnell & Price| ( |1995| ) and 
|Price & Bate|(|2007| > and parallelised using both OpenMP and MPI 
(see |Meru"& Bate 201 lb for details). Our simulations with a grid- 
based code are carried out using the Fast Advection in Rotating 
Gaseous Objects (FARGO) two-dimensional fixed polar hydrody- 
namics code (Masset 2000; Baruteau & Masset 2008a b). 

We include the heating effects in the disc due to work done on 
the gas and artificial viscosity. The cooling in the disc is taken into 
account using the cooling parameter, p (equation[2](, first proposed 
by |Ga mmie (2001) which cools the gas on a timescale given by 
equation [3] For the SPH simulations carried out in this paper, we 
ensure that the timestepping is limited by the following timestep 
criterion (in addition to the Courant condition, the force condition 
and the viscous timestep condition; see Monaghan 1992 ): 

A*<og, (5) 

where C = 0.3. Meru & Bate (2011b) show that this condition is 
adequate to ensure that the fragmentation results are not affected 
by the timestepping imposed. For the simulations performed using 
FARGO, the timestep constraint (using C = 1) is also included 
for all simulations involving p < 6. This constraint appears as an 
additional term in the denominator of equation 15 of Masset (2000). 
We have verified that for larger values of p the effect of including 
this timestepping constraint is negligible. 

To model the shocks in the discs, both codes use artificial vis- 
cosity. The SPH code uses the artificial viscosity method described 
by |Chow & Monaghan| (1997} and |Price & Monaghan| |20 04], 
the implementation of which is summarised in equations |A1| |A2| 
and |A10| (see Appendix |A| for details). The artificial viscosity is 
controlled by the parameters asPH and /?sph- FARGO uses the 
|von Neumann & Richtmyer ( 1950) artificial viscous pressure with 
parameter q. We use the default values for the SPH and FARGO 
artificial viscosity parameters of (q;sph, /^sph) = (0.1,0.2) and 
q — 1.41, respectively. Where specified, we also vary the amount 
of SPH and FARGO artificial viscosities to investigate their effects 
on the fragmentation boundary. 

FARGO's grid is set up to use linear spacing in the azimuthal 
direction and logarithmic spacing in the radial direction. We use 
open boundary conditions at both the inner and outer radial bound- 
aries and use a fixed gravitational softening length of 3 x W~ 4 H, 
where H is the vertical scaleheight of the disc, in all the FARGO 
simulations. 
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3 THE EFFECT OF NUMERICAL VISCOSITY 



In hydrodynamics codes, artificial viscosity is frequently applied 
to correctly capture shocks and to avoid post-shock oscillations. In 
trying to understand the evolution of gravitationally unstable discs, 
equation[4]seeks to link the dissipation rate in the disc to the cooling 
rate of the disc. However, equation [4] is derived for a steady-state 
disc which assumes that gravitational instability is the only heating 
source for the disc. In reality, in numerical simulations, there will 
be additional heating due to numerical dissipation. 

In SPH, the artificial viscosity typically includes both a linear 
term controlled by osph and a quadratic term controlled by /3sph- 
The linear term provides a bulk viscosity which dissipates kinetic 
energy as particles approach each other to reduce particle oscilla- 
tions following a shock, while the quadratic term is primarily im- 
portant to stop particle interpenetration. FARGO, on the other hand, 
only contains a quadratic term controlled by the artificial viscosity 
parameter, q, and provides a bulk viscosity. 

The dissipation from the bulk viscosity in shocks (such as 
those generated by gravitational instability) is physical. In a grav- 
itationally unstable disc, this provides the «gi term. However, in 
these discs the artificial viscosities will also provide some shear 
viscosity, the heating effects of which are not accounted for in equa- 
tiong] 



3.1 SPH artificial viscosity 



Monaghan (1985) showed that in the continuum limit, the qsph 
component mimics a Navier-Stokes viscosity with bulk and shear 
coefficients proportional to the resolution length (see Meglicki et al. 
1993 and Appendix |A]>. This has been confirmed numerically (e.g. 



Artymowicz & Lubow|1994}|Lodato & Price|2010|l. The shear vis- 



cosity contributions from the linear and quadratic SPH terms due 
to the artificial viscosity can also be compared to the Sha kura &| 
Syunyaev ( 1973 1 viscosity of the form v = assc B H. For the SPH 
calculations discussed in this paper, where the viscosity is only ap- 
plied between approaching particles, assuming a Keplerian disc it 
can be shown that (see Appendix |A2[|B2| and|C]l 



Q!SS,lin 



31 h 
525 QSPH tf 



(6) 



9 / 31 h 9 a 



(8) 



Note that while the dissipation due to the quadratic term is often ig- 
nored when comparing the viscosity in SPH simulations of a-discs, 
we show in Sections |5.3.5| and Appendix [C] that its contribution 
can be non-negligible. In a numerical simulation, it is this heating 
rate that must be balanced by the imposed cooling (equation[3} for 
the disc to settle into a quasi-steady state. Thus, equation [4] can be 
rewritten as 



9 -y(7 — 1) (a GI + a S S,lin + OSS, quad) 



(9) 



For the disc to fragment, the combined heating must be insufficient 
to balance the cooling so that 



/3crit — 



1 



' 7(7 — 1) ( a GI,crit + QSS,lin,crit + asS.quad.crit ) ' 

(10) 

where aGi.crit is the true value of the gravitational stress that a disc 
can support before fragmenting and a S s,iin, C rit and a S s, q uad,crit 
are the contributions to the heating from the artificial viscosity that 
allows a disc to fragment once /3 < /3 cl -i t for any one particular 
resolution. 

Now, for a particular cooling self-gravitating disc calculation, 
let us suppose that there is some maximum gravitational stress that 
can be produced by the disc beyond which it will fragment. In this 
case, a G i,crit will be a constant, but Qss,iin,crit and a S s, q uad,crit 
will decrease with increasing resolution. If ass.iin and ass, quad 
obey equations[6]and[7] then for a set of simulations with increasing 
resolution 



/3crit — ^ 



9 7(7 



31 h , 9 a ( h 

a GI , crit +^a SP H^+C^SPH 

(11) 



and 



where rj and £ are constants which we expect to equal 1 . 

In order to compare equation [TT] with the results of SPH sim- 
ulations, we need to determine h/H just before the fragmentation 
sets in. Assuming that the disc fragments when Q ~ 1, using equa- 
tion [T] and noting H = c B /Q, and Q 2 — GM*/i? 3 , we have 



OSS, quad 



70tt 



/SSPH 



(7) 



where ass.iin and ass, quad are the contributions from the linear 
and quadratic terms, respectively, and h is the particle smoothing 
length. Note that the coefficients in this SPH code are marginally 
different to some other SPH implementations. In SPH codes em- 
ploying the older Monaghan & Gingold ( 1983 ) formalism, the co- 
efficients would be 1/20 and 3/(357r) for the linear and quadratic 
terms, respectively (see Appendix | A 1 1 |B 1 1 and |C| f or details). 

For a Shakura & Syunyaev ( 1973 1 disc model, the dissipation 
rate per unit mass is given by |assc a Q. In a purely gravitationally 
unstable disc the dissipation rate may be parametrized |aGiCsfi. 
However, using SPH we expect an additional heating due to nu- 
merical dissipation given by | (ass.iin + ass,quad)cf fi. Thus, the 
combined heating rate per unit mass is expected to be 



H 



(12) 



where R is the radius in the disc and M* is the mass of the central 
object. The smoothing length in an SPH simulation for a disc that 
is resolved (i.e. h < H, which is true for all the simulations pre- 
sented here a t the radius of fragmentation) can be estimated, using 
equation 3 of Price & Bate ( 2007 1 with density, p « =^7 , as 



2Hmr 



1/3 



(13) 



where m p is the mass of an SPH particle. We use constant mass 
SPH particles and so the mass of an SPH particle is the disc mass 
divided by the number of SPH particles, m p = Md/N pa , Tt . The 
ratio of the smoothing length to disc scaleheight can then be ap- 
proximated to be 
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h 

H 



1.2 



2M 2 M d 
n 2 N, 



part 



1/3 



(14) 



In Section [5. 3.4| we verify that this is indeed the case for steady- 
state marginally stable discs that have a Toomre parameter, Q w 1. 
We therefore expect that 



Pcrit — 



9 7(7-1) 



31 1.2 /2M 2 M d \ 1/3 

CKGI.crit + Vf^ I ~^1Tt ] «SPH 



'525 E.R 2 \n 2 N p 



u 9 (1-2) 2 A gAgMd V 

+ S07T E 2 i? 4 U 2 ^part7 PSPH 



(15) 



i.e. we have three unknowns: aGi,crit, f] and £, since we can de- 
termine /3 cr it for any one resolution. If 77 and £ are unity, then the 
second and third terms in the denominator in equationfT31each have 
values of O(10 -3 ) for the discs studied by |Rice et al.| {2005} and 
|Meru & Bate]j201 la{ l simulated with 250,000 particles (using the 
parameters described in Section [4] and given that the fragmenta- 
tion occurs in the outer parts of the disc - see Meru & Bate 201 lb| >. 
The contribution from these terms are approximately a factor of 
O(10) smaller than the original estimate of «Gi,crit ~ 0.06 (Gam- 
mie ptjQT] |Rice et al.|2005> . Thus it was assumed in earlier SPH 
studies that the heating due to artificial viscosity would be negli- 
gible compared to the dissipation due to gravitational instabilities 
(see Appendix A of Lodato & Rice 2004 1. Note, however, that if the 
SPH artificial viscosity plays a significant role then the qsph term 
scales linearly with the smoothing length, h, such that the conver- 
gence of /3 C iit towards the true value is expected to be first order in 
h (i.e. very slow as the numerical resolution is increased). On the 
other hand, if the dominant term is the /3sph term, the convergence 
will be faster since the SPH artificial viscosity scales quadratically 
with the smoothing length. 



3.2 fargo artificial viscosity 

Most grid-based hydrodynamical codes are second order and do not 
have a linear viscosity. Therefore, one would expect that their rate 
of convergence towards the true value of /3 cr it will be second order 
in spatial resolution and thus possibly faster than SPH codes. In- 
deed, FARGO uses the |von Neumann & Richtmyer|([1950|> art ificial 
viscous pressure given by (see also |Bodenheimer et al.|2007) 



Using equation [T7] and equating equation [16] to the shearing force 
per unit area as defined in equation |B 1 j yields a kinematic viscosity 
due to the artificial viscosity given by 



^av,fargo = Q (Ax) 



R 



dQ. 
dR 



(18) 



where the final approximation assumes a Keplerian flow. This 
gives a |Shakura & Syunyaev| {1973| > type viscosity of the form 
v = assCsH with 



QSS, fargo = »Q 



Ax 

If 



(19) 



Analogous to the derivation of equation [TT] for SPH artificial 
cosity this yields 



97(7-1) 



1 

CKGI.crit + £§ q 2 



(20) 



where we expect that £ is unity. Substituting for the scaleheight 
using equation [T2] yields a formula for /3 C rit that is equivalent to 
equation[L5]for a grid-based code: 



r 4 1/ 3 2 {AxM. 

^ = 9 7^1) i QGI ' crit+C 2 9 {^W 



(21) 



Since the FARGO simulations use a logarithmic grid we take the 
cell size at a radius of R — 22 au, i.e. close to the edge of the disc 
where we would expect fragmentation to occur when the cooling 
timescale is close to the critical value for any one particular resolu- 
tion. We find that the cell size at this radius scales as 



Ax ! 



125 

nL. 



(22) 



where N cc \\ s is the total number of cells used in the simulation. We 
therefore have two unknowns: QGi.crit and £, since we can deter- 
mine /3 C rit at any one resolution. 

Equation|20|shows that if the artificial viscosity plays a signif- 
icant role in the dissipation in the disc, the convergence is expected 
to be second order in spatial resolution, i.e. potentially faster than 
with SPH. 



<? 2 p(Ax) 2 



dv 
dx 



(16) 



where Ax is the cell size, p is the density and q = I /Ax is a 
constant which indicates the number of grid cells over which the 
shock is spread and whose value is dependent on the numerical 
scheme and is usually 0.05 q ^ 2 (I indicates the strength of the 
artificial viscosity). This is a bulk viscosity. In a cylindrical code, if 
the gas travels in circles, there should be no shear viscosity at all. 
However, in a gravitationally unstable disc, this will not be the case 
and there will be some shear viscosity (that arises from the bulk 
viscosity) which we expect will scale in roughly the same way, i.e. 
proportional to the square of the size of the grid cell. Assuming the 
shear rate, \dv/dx\, is approximately Keplerian, then 



dv 

dx 



R 



dR 



(17) 



3.3 Artificial viscosity effects with resolution 

In an SPH code when the number of particles increases, h/H 
is reduced and therefore as the resolution increases, ass,iin and 
oss, quad, and thus the heating due to artificial viscosity, tend to 
zero. Similarly, as the resolution is increased in a grid-based code, 
the cell size decreases for any one problem, and the contribution to 
the dissipation from the artificial viscosity decreases. In the limit 
of infinite resolution, equations [TT] and [20] will return to equation 
[4] But for a finite resolution, the value of /? C rit obtained from a nu- 
merical simulation should always be smaller than the true value. 

It is also important to note that simply reducing the value of 
Qsph, /3sph or q is not necessarily a sufficient way in which to de- 
crease the numerical dissipation and obtain the true value of /3 cr it- 
By inspection of equations [TT] and [20] we might naively assume 
this to be the case. However, reducing these values may mean that 
shocks are treated inaccurately. For example, the shocks may not 
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be spread over a large enough lengthscale to model them numeri- 
cally and/or there may be post-shock oscillations that are eventually 
damped, resulting in dissipation. Lodato & Price (2010) showed us- 
ing SPH that setting qsph = counter-intuitively led to a larger 
amount of dissipation. |Price & F ederrath ( 2010) found that if an 
adequate value of /3sph is not used, particle interpenetration may 
occur. In their case, they stated that this makes very little difference 
to their dissipation rate since their linear term dominates almost ev- 
erywhere. However, their simulations explored a different regime 
to that being explored here and with a different artificial viscos- 
ity scheme. Their simulations employed the Morris & Monaghan 
( 1996 ) artificial viscosity switch where the value of «sph ranges 
between 0.05 and 1.0 (the higher value being implemented close 
to shocks). Their simulations were of high Mach number shocks 
(./£ = 10) and so a large part of their simulations would require 
the use of the higher value of qsph and thus this would dominate 
the dissipation. The simulations performed by |Rice et al.|p005) 
and Meru & Bate (201 la) used a fixed value of «sph = 0.1 - such 
a low value opens up the possibility of the quadratic term being im- 
portant and therefore decreasing the value of /?sph may affect the 
overall dissipation rate and thus the fragmentation outcome. 

In FARGO, if q is set to zero, there will be no controlled numer- 
ical dissipation (e.g. to capture shocks or other disturbances at the 
grid scale). However, the code will still have some level of numer- 
ical diffusion and dissipation which is not controllable other than 
that it too should decrease with increasing resolution. 



4 SIMULATIONS 

The disc and star properties used to carry out the simulations in this 
paper are exactly the same as those used by Rice et al. (2005) and 
|Meru&Bate|j2011a| l: a O.1M disc surrounding a 1M Q star. The 
SPH simulations span a radial range, 0.25 ^ R ^ 25 au, while 
the FARGO simulations span a radial range, 1 ^ R ^ 25 au (only 
marginally different to the SPH simulations for numerical reasons). 

The initial surface mass density and temperature profiles are 
E oc and T oc R~ x ^ 2 , respectively, and the temperature is 
normalised so that the minimum initial Toomre stability value at 
the outer edge of the disc, Q m i n = 2. The discs are modelled with 
a ratio of specific heats, 7 = 5/3. 

Table [T] shows a summary of the initial SPH simulations and 
the key fragmenting results carried out by Meru & Bate ( 2011a) 
(obtained from their Table 1) as well as those in this paper (bold 
text). We supplement the Meru & Bate (201 la) results by carrying 
out an additional simulation using 2 million particles with /3 = 
9 and three additional simulations using 16 million particles with 
P = 12, 15 and 20. 

Table [2] summarises the initial simulations carried out using 
FARGO and the key fragmentation results. We perform simulations, 
using q — 1.41, at five different resolutions and determine the crit- 
ical value of p at each of these resolutions. The lowest resolution 
simulations are carried out using 768 and 256 grid cells in the az- 
imuthal and radial directions, respectively. We then increase the 
linear resolution by factors of 2, 4, 8 and 16 in both the azimuthal 
and radial directions. 

The simulations were run either for at least 6 outer rotation pe- 
riods (ORPs) or until the discs fragmented. Fragments are defined 
as regions whose surface mass densities are at least two orders of 
magnitude denser than their surroundings. In addition, we ensure 
that the fragments survive for at least one rotation to verify that 
they do not shear apart. 



Simulation name 


No of particles 




Fragmented? 


31k-beta2 


31,250 


2.0 


Yes 


31k-beta2.5 


31,250 


2.5 


Yes 


31k-beta3 


31,250 


3.0 


Yes 


31k-beta3.5 


31,250 


3.5 


No 


31k-beta4 


31,250 


4.0 


No 


250k-beta5 


250,000 


5.0 


Yes 


250k-beta5.5 


250,000 


5.5 


Yes 


250k-beta5.6 


250,000 


5.6 


No 


250k-beta6 


250,000 


6.0 


No 


250k-beta6.5 


250,000 


6.5 


No 


250k-beta7 


250,000 


7.0 


No 


250k-beta7.5 


250,000 


7.5 


No 


2m-beta5.5 


2 million 


5.5 


Yes 


2m-beta6 


2 million 


6.0 


Yes 


2m-beta6.5 


2 million 


6.5 


Yes 


2m-beta7 


2 million 


7.0 


Yes 


2m-beta8 


2 million 


8.0 


Yes 


2m-beta9 


2 million 


9.0 


No 


2m-betal0 


2 million 


10.0 


No 


2m-betal0.5 


2 million 


10.5 


No 


2m-betal 1 


2 million 


11.0 


No 


2m-betal5 


2 million 


15.0 


No 


16m-betal0 


16 million 


10.0 


Yes 


16m-betal2 


16 million 


12.0 


No 


16m-betal5 


16 million 


15.0 


No 


16m-betal8 


16 million 


18.0 


No 


16m-betal0 


16 million 


20.0 


No 



Table 1. Table showing the SPH simulations carried out by Meru & Bate| 
(2011a) , as well as the supplementary SPH simulations performed in this 
paper (bold text), and the key fragmenting results. The simulations are per- 
formed using («sph, /3sph) = (0.1,0.2) 



To investigate the effects of the different components of ar- 
tificial viscosity in SPH, we carry out a number of simulations 
where we vary the values of qsph and Psph (see Table [3] for de- 
tails). We perform a suite of simulations using 250,000 particles. 
Firstly, we set qsph = 0.1 and vary the value of Psph (Table[3] 
top section). Secondly, we set Psph ~ 2.0 and vary the value of 
Q!Sph (Table[3] middle section). We then carry out simulations with 
(o!Sph,,3sph) = (0.1,2.0) using 31,250, 250,000 and 2 million 
particles to determine the effect that changing the value of Psph 
has on the fragmentation boundary at each of these resolutions (Ta- 
ble^ bottom section). 

In addition, we perform a number of SPH simulations without 
self-gravity (see Table|4]for details) using various values of «sph, 
Psph to compare the measured dissipation to the analytically ex- 
pected values in equation[8](also see Appendix [C]). 

Finally, we investigate the effects that artificial viscosity in 
FARGO has on the critical cooling timescale by varying the value 
of q in equation[T6]between and 2.5 (see Table[5j for discs mod- 
elled using 786,432 grid cells (512 and 1536 cells in the radial and 
azimuthal directions, respectively). We then carry out simulations 
using an artificial viscosity parameter, q = 0.5, at all but the lowest 
resolutions considered in this paper to determine the effect that this 
has on the fragmentation boundary. 
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Figure 2. Surface mass density rendered image of two identical simulations carried out using fi = 10 (left panel) and /3 = 12 (right panel), using 16m SPH 
particles. Fragmentation occurs for /3 = 10 but not for /3 = 12. These simulations are performed with (crspH, /9sph) = (0.1, 0.2). 
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Figure 1. Graph of against resolution of the non-fragmenting (open 
squares) and fragmenting (solid triangles) SPH simulations. This figure 
contains the results presented by Meru & Bate (2011a their Figure 3) as 
well as the new simulations highlighted in TablefT] The solid line, obtained 
by fitting equation [23] shows a dividing line between the fragmenting and 
non-fragmenting cases and the grey region is where fragmentation can take 
place. The graph shows clear evidence of convergence of results with in- 
creased resolution. These simulations are carried out with (qsph, /?SPh) 
= (0.1, 0.2). The convergence rate is first order with spatial resolution. The 
dotted line (which coincides well with the solid line) is obtained by fitting 
equation|15| 



5 RESULTS 

5.1 The convergence rate of /3 cr it with SPH 

Table[T]summarises the results of the SPH simulations, using aspH, 
/3 SPH ) = (0.1, 0.2), carrie d out by|Meru & Bate| < |2011a| and those 
performed for this paper. [Meru & Bate (201 la i found borderline 
simulations which they defined to be discs which showed signs 
of fragmentation but the fragments sheared apart rapidly (within 
1 ORP) and no further signs of fragmentation were seen. We find 
that borderline simulations can in fact range a span of fi values. 
However, since ultimately they are discs that do not end up frag- 
menting, we now simplify this terminology and refer to them as 
non-fragmenting simulations. Figure [T] shows a summary of the 
SPH results. This figure is the same as Figure 3 of |Meru & Bate| 
l |2011a| > but with the the addition of the new SPH results presented 
in this paper. |Meru & Bate| ( |2011a| > found no evidence for con- 
vergence, but the addition of the new high resolution calculations 
now provides evidence for a very slow convergence of /3 cr it with 
increasing resolution. Figure|2]shows two of the highest resolution 
SPH simulations (16 million particles) carried out using fi = 10 
and fi = 12. It can clearly be seen that at this resolution, fragmen- 
tation occurs for fi = 10 but not for fi = 12. To estimate the rate 
of convergence we fit a formula of the form 



(23) 



l + A/ ff 

where I is the linear spatial resolution, A is a constant and a is the 

_ i 

convergence rate. For SPH, we simply take I oc N pg ^ t . We fit this 
formula to the values of fi in Figure [T] that lie half way between 
the lowest non-fragmenting value of fi and the highest fragment- 
ing value of fi for each numerical resolution, i.e. the fragmentation 
boundary. We find that a good fit is obtained with/3 cr i t = 15.6±1.0 
and a = 1.08 ± 0.05. The value of a shows that the rate of con- 
vergence is first order in spatial resolution. For the benefit of un- 
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Simulation name 


No of 
radial cells 


No of 

azimuthal cells 


P 


Frag 


107k rplk-hptaO S 

17 / A._CClla ~UCLa\J.J 


256 


768 


0.5 


Yes 


1 Q71r <-plk hptal 
17 /K_ccll&-UcLai 


-JO 


768 


j 


Yes 


1 07V rplk-hpfa? 


256 


768 


2 


No 


1 Q7V /-plk hptn^ 
IV /K-CcllS-UcLdj 


256 


768 




No 


1 Q7V /-plk-hptaA 

17 / A._CClla UCLtV-r 


256 


768 


4 


No 


7SAV /-plk hpta^ 


SI 7 




2 


Yes 


7R6V rplk-hpfa^ S 

/ OUd._LCllft UCLtlJ .J 


5 12 


1536 


3 5 


Yes 


7R6V <-plk-hpfa4 5 

/ OU(i._CC11l> UCLtlT^ J 


512 


1536 


4.5 


Yes 


7R6k rplk-hpia5 


512 


1536 


g 


Yes 


/-plk hptaS S 


S 1 9 
J 1 z 




g g 


Yes 


78fM- /-pile hpla^ 
/ OOK_Ccll„-UcLdO 


S 1 7 


1 


u 


No 


IRfciV /-plk hpfalf) 
/ OUK_Ccllft UcLdlU 


S 1 7 


1 

1 J JO 


i n 


No 


3. lm_cells-betalO 


1024 


3072 


i n 


Yes 


3. 1 m_cells-betal 2 


1024 


3072 


1 2 


Yes 


a 1 m /-p 1 k _ nplci 1 a 
J. 1 111_CC11„ UCld.1 J 


1024 


3072 


1 3 


No 


1 3m_cells-beta 1 1 


2048 


6144 


1 1 


Yes 


1 Sm /— plk«r\ptn 1 ZL 
1 J111_CC11_ UCld L'-r 


2048 


6144 


14 


Yes 


1 3m_cells-beta 1 5 


TO/1 9 


A 1 A A 
1*1**' 




Vac 

I es 


13m_cells-betal6 


2048 


6144 


16 


Yes 


13m_cells-betal8 


2048 


6144 


18 


No 


50m_cells-betal8 


4096 


12288 


18 


Yes 


50m_cells-beta20 


4096 


12288 


20 


Yes 


5 0m_cells-beta22 


4096 


12288 


22 


Yes 


50m_cells-beta24 


4096 


12288 


24 


No 



Table 2. Table showing the simulations carried out using FARGO and the 
key fragmenting results. The simulations are performed using the artificial 
viscosity parameter, q = 1.41. 

derstanding the results presented in Section |5.3.3| (which only have 
data points at the lowest three resolutions) we fit equation[23]to the 
data presented in Figure^but exclude the highest resolution simu- 
lations. We find /3 cr i t = 17.4 and a = 1.03, so excluding the last 
point does not alter the fit significantly due to the slow convergence 
rate of SPH. 

In addition, we also fit equation [15] to this data. We find that 
aGi.cfit = 0.024 ± 0.001, rj = 21.1 ± 1.3 and £ = 1.7 ± 1.5. 
In the limit of infinite resolution, this value of the critical gravita- 
tional stress obtained is equivalent to a critical cooling timescale 
of /3 C rit ~ 17. However, we point out that while £ is reasonably 
close to unity, the value of rj is very large. This is suggestive of an 
additional source of dissipation present in the simulations over and 
above what we expect from artificial viscosity in a shear dominated 
disc. 

5.2 The convergence rate of /? cr it with fargo 

Table [2] and Figure [3] summarise the results using the grid-based 
code, FARGO. As with the SPH results, we also see that as the res- 
olution increases numerical convergence does appear to take place. 
Figure [4] shows the surface mass density rendered images of the 
discs modelled using FARGO at resolutions of 786,432, 3.1 million, 
13 million and 50 million grid cells. The top two panels show that 
using a cooling timescale of j3 = 10, the discs do not fragment 
using 786,432 grid cells but when the resolution is increased to 3. 1 
million grid cells, fragmentation is seen. Similarly, using a cooling 
timescale of /3 = 18, the disc modelled with 13 million grid cells 
does not fragment whereas that modelled at the higher resolution 
of 50 million grid cells does fragment. 

Again, we use equation|23]to estimate the rate of convergence. 
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Figure 3. Graph of /? against resolution of the non-fragmenting (open 
squares) and fragmenting (solid triangles) FARGO simulations carried out 
using q = 1.41. The solid line, obtained by fitting equation |23| shows a 
dividing line between the fragmenting and non-fragmenting cases and the 
grey region is where fragmentation can take place. The graph shows clear 
evidence of convergence of results with increased resolution. The conver- 
gence rate is second order with spatial resolution. The dotted line (which 
coincides well with the solid line) is obtained by fitting equation [2T] 

_ i 

For FARGO, we simply take I oc N cc ^ s , where 7V cc n s is the number 
of grid cells. Note that the linear resolution is inversely proportional 
to the square root of the number of cells because the calculation is 
two dimensional. We fit this formula to the fragmentation boundary 
in Figure [3]( as done for the SPH results. We find that a good fit 
is obtained with /? crit = 22.3 ± 2.3 and a = 2.03 ± 0.36. The 
value of a shows that the rate of convergence is second order in 
spatial resolution. We then fit the data using equation |2T] and find 
that aci.crit = 0.018 ± 0.001 and £ = 0.87 ± 0.08. In the limit of 
infinite resolution, this value of aGi.crit is equivalent to a critical 
cooling timescale, /3 C rit ~ 22 (using equation H), similar to the 
value obtained using equation [23] 

Thus the value of /3 cr jt converges more rapidly using FARGO 
than SPH. In Section [3~2] we note that if artificial viscosity plays a 
significant role in the determination of /3 cr it then FARGO might be 
expected to display second-order convergence since it only applies 
a quadratic artificial viscosity. On the other hand SPH includes both 
linear and quadratic artificial viscosities. If the linear term is domi- 
nant, this may lead to first-order convergence. This implies that arti- 
ficial viscosity may be significant in determining /3 cr it . We note that 
FARGO appears to converge to a higher value of /3 cr i t than SPH, but 
this result may also be caused by the different artificial viscosities. 
Therefore, in the following sections we investigate the dependence 
of /3 C rit on the artificial viscosities applied in both codes. 

5.3 The effect of SPH artificial viscosity on convergence 

In Section[3]we present analytical arguments that suggest that arti- 
ficial viscosity may play a role in the numerically determined value 
of the critical cooling timescale. We show that the contribution to 
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Figure 4. Surface mass density rendered images of four simulations carried out using FARGO with cooling timescales of @ = 10 (top panel) and /3 = 18 
(bottom panel) using q = 1.41. For the simulation with 786,432 grid cells (upper left panel), fragmentation does not occur when modelled with = 1Q but 
fragmentation is seen when the resolution is increased to 3.1 million grid cells (upper right panel). Similarly, when the resolution is increased further to 13 
million grid cells (lower left panel), fragmentation is not seen using a cooling timescale of /3 = 18 whereas when the same simulation is carried out using 50 
million grid cells (lower right panel) fragmentation is seen. 



the dissipation due to the artificial viscosity is expected to decrease 
with increasing resolution (Appendix[C]and equation[8j. Therefore, 
if the slow convergence can be attributed to SPH artificial viscosity, 
this may be the reason why Meru & Bate (201 la) found that /3 cr it 
increases with increasing resolution and is a plausible explanation 
as to why the results presented in Section |5~T| show a slow conver- 
gence. We test the role that artificial viscosity plays on the frag- 
mentation of self-gravitating discs by varying the values of osph 
and /3sph separately. 



5.3.1 The effect o//?sph on the critical cooling timescale 

In Section [5~T| we show that the first-order convergence seen for the 
SPH results suggests that the «sph term may be responsible. How- 
ever, as mentioned in Section [33] if the optimum value of /3sph 
is not used (i.e. a value that minimises numerical dissipation) addi- 
tional dissipation may occur and affect the fragmentation boundary. 
Therefore, while not immediately obvious from the results in Sec- 
tion |5T| the value of the /3sph term may affect the fragmentation 
conclusions. Table[5](top panel) and Figure[3]summarise the results 
of the simulations carried out to investigate what effect the value 
of Psph has on the critical cooling timescale using 250,000 parti- 



cles and maintaining a fixed value of osph = 0.1. It can be seen 
that the shape of the fragmenting/non-fragmenting boundary line 
appears to follow a somewhat S-shaped curve. At high values of 
Psph, any potential particle interpenetration is appropriately dealt 
with as changing the value of /?sph from 2 to 4 has no effect on the 
critical cooling timescale. As /?sph is reduced, particle interpen- 
etration and additional particle velocity dispersion can occur since 
the appropriate amount of the quadratic term of the artificial viscos- 
ity is not used. Eventually, these oscillations are damped down by 
the q:sph term resulting in dissipation. Consequently, a more rapid 
cooling is required to overcome the additional dissipation result- 
ing in smaller critical cooling values. At very low values of /3sph 
the critical cooling timescale remains the same. This may be due 
to one of two reasons: 1) a "saturation" of additional oscillations 
occurs such that reducing /?sph further does not increase the dissi- 
pation - by this we mean that the cause of the oscillation (i.e. the 
incorrect modelling at the edge of the shock front and particle in- 
terpenetration) is so high that any reduction in /3sph cannot cause 
more oscillation to occur; or 2) at such low values of /3sph, the 
linear artificial viscosity term dominates the dissipation such that 
any additional particle interpenetration does not increase the over- 
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Simulation name 


No of particles 


^SPH 


/^SPH 


P 


Fragmented? 


250k-betaSPH0.1-beta5 


250,000 


0.1 


0.1 


5.0 


Yes 


250k-betaSPH0.1-beta5.5 


250,000 


0.1 


0.1 


5.5 


Yes 


250k-betaSPH0.1-beta5.6 


250,000 


0.1 


0.1 


5.6 


No 


250k-betaSPH0.1-beta6 


250,000 


0.1 


0.1 


6.0 


No 


250k-beta5 


250,000 


0.1 


0.2 


5 


Yes 


250k-beta5.5 


250,000 


0.1 


0.2 


5.5 


Yes 


250k-beta5.6 


250,000 


0.1 


0.2 


5.6 


No 


250k-beta6 


250,000 


0.1 


0.2 


6.0 


No 


250k-beta6.5 


250,000 


0.1 


0.2 


6.5 


No 


250k-beta7 


250,000 


0.1 


0.2 


7.0 


No 


250k-beta7.5 


250,000 


0.1 


0.2 


7.5 


No 


250k-betaSPH0.4-beta6 


250,000 


0.1 


0.4 


6.0 


Yes 


250k-betaSPH0.4-beta6.5 


250,000 


0.1 


0.4 


6.5 


Yes 


250k-betaSPH0.4-beta6.8 


250,000 


0.1 


0.4 


6.8 


No 


250k-betaSPH0.4-beta7 


250,000 


0.1 


0.4 


7.0 


No 


250k-betaSPHl-beta6.5 


250,000 


0.1 


1 


6.5 


Yes 


250k-betaSPHl-beta6.8 


250,000 


0.1 


1 


6.8 


Yes 


250k-betaSPHl-beta7 


250,000 


0.1 


1 


7.0 


No 


250k-betaSPH2-beta4 


250,000 


0.1 


2 


4 


Yes 


250k-betaSPH2-beta5 


250,000 


0.1 


2 


5 


Yes 


250k-betaSPH2-beta6 


250,000 


0.1 


2 


6 


Yes 


250k-betaSPH2-beta7 


250,000 


0.1 


2 


7 


Yes 


250k-betaSPH2-beta8 


250,000 


0.1 


2 


8 


Yes 


250k-hetaSPH2-hpta8 5 


250 000 


0.1 


2 


8.5 


No 


250k-hetaSPH2-hetn9 


250 000 


0.1 


2 


9 


No 


250k-hetaSPH2-hetn 1 


250 000 


0.1 


2 


10 


No 


Z J UK- Uc Lao rlH- - DC Lao 


?sn nnn 

jLJKJ,\AJU 


n 1 

U. 1 


4 


R n 


ICS 


Z JUK-UcLaor^ri i t-Dfc:ldO. J 


?sn non 


n 1 

V. I 


4 


O.J 


1NO 


250k-alphaSPH0.05-beta7 


250,000 


0.05 


2 


7 


Yes 


250k-alphaSPH0.05-beta8 


250,000 


0.05 


2 


8 


No 


250k-alphaSPH().05-beta9 


250,000 


0.05 


2 


9 


No 


250k-betaSPH2-beta4 


250,000 


0.1 


2 


4 


Yes 


250k-betaSPH2-beta5 


250,000 


0.1 


2 


5 


Yes 


250k-betaSPH2-beta6 


250,000 


0.1 


2 


6 


Yes 


250k-betaSPH2-beta7 


250,000 


0.1 


2 


7 


Yes 


250k-betaSPH2-beta8 


250,000 


0.1 


2 


8 


Yes 


250k-betaSPH2-beta8.5 


250,000 


0.1 


2 


8.5 


No 


250k-betaSPH2-beta9 


250,000 


0.1 


2 


9 


No 


250k-betaSPH2-betalO 


250,000 


0.1 


2 


10 


No 


250k-alphaSPH0.2-beta7 


250,000 


0.2 


2 


7 


Yes 


250k-alphaSPH0.2-beta7.5 


250,000 


0.2 


2 


7.5 


Yes 


250k-alphaSPH0.2-beta8 


250,000 


0.2 


2 


8 


No 


250k-alphaSPH0.5-beta6 


250,000 


0.5 


2 


6 


Yes 


250k-alphaSPH().5-beta6.5 


250,000 


0.5 


2 


6.5 


Yes 


250k-alphaSPH().5-beta7 


250,000 


0.5 


2 


7 


Yes 


Z JUK-dipildi3J^rlW. J-Ucld / .J 


9sn non 

ZJU, \J\JU 


n s 

U.J 


2 


7 S 
/ .J 


1NO 


Z JUK-dipildk3I^n 1-UcLdJ 


9sn nnn 


I 


2 


J 


ies 


Z JUK-dipildk3X^rl 1-UcLdO 


9sn nnn 


1 


2 


f. 
o 


Vps 

ies 


Z JUK-dipildi3J^rl 1-UcLdO. J 


9sn nnn 


] 


2 


6 S 

O.J 


Ypc 

ies 


Z JUK-dipildk3J^rl 1-UcLd / 


9sn nnn 


] 


2 


7 
1 


1NO 


31k-betaSPH2-beta3.5 


31,250 


0.1 


2 


3.5 


Yes 


31k-betaSPH2-beta4 


31,250 


0.1 


2 


4 


Yes 


31k-betaSPH2-beta4.5 


31,250 


0.1 


2 


4.5 


Yes 


31k-betaSPH2-beta5 


31,250 


0.1 


2 


5 


No 


250k-betaSPH2-beta4 


250,000 


0.1 


2 


4 


Yes 


250k-betaSPH2-beta5 


250,000 


0.1 


2 


5 


Yes 


250k-betaSPH2-beta6 


250,000 


0.1 


2 


6 


Yes 


250k-betaSPH2-beta7 


250,000 


0.1 


2 


7 


Yes 


250k-betaSPH2-beta8 


250,000 


0.1 


2 


8 


Yes 


250k-betaSPH2-beta8.5 


250,000 


0.1 


2 


8.5 


No 


250k-betaSPH2-beta9 


250,000 


0.1 


2 


9 


No 


250k-betaSPH2-betalO 


250,000 


0.1 


2 


10 


No 


2m-betaSPH2-betal2 


2 million 


0.1 


2 


12 


Yes 


2m-betaSPH2-betal4 


2 million 


0.1 


2 


14 


No 



Table 3. Table showing the simulations carried out to investigate how the fragmentation boundary changes with the SPH artificial viscosity parameters (i) 
^^'QJ^^aMMr^^^u^rng 250,000 particles (upper panel), (ii) /3 S ph = 2.0 and varying a S pn using 250,000 particles (middle panel) and (iii) 
(agpH. /3sph) = (0-1. 2.0) at different resolutions (lower panel). The key fragmenting results are also indicated. 
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Figure 5. Graph of f3 against /3sph of the non-fragmenting (open squares) 
and fragmenting (solid triangles) simulations carried out using 250,000 par- 
ticles and aspH = 0.1. The solid line, included by eye, shows a dividing 
line between the fragmenting and non-fragmenting cases and the grey re- 
gion is where fragmentation can take place. The graph shows an S-shaped 
curve: any additional particle oscillation that may exist is stopped using 
i^SPH ~ 2 since the effect on the critical cooling timescale does not change 
above this value; at lower values of /3gpH, the critical cooling timescale is 
smaller as the additional particle oscillation results in excess dissipation that 
needs to be overcome before fragmentation can take place. At very low val- 
ues of /3sph either so much particle oscillation occurs (either at the edge 
of the shock front or due to particle interpenetration) or that the osph term 
dominates, that the effects of lowering /?sph does not result in more dissi- 
pation. 
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Figure 6. Graph of /3 against osph of the non-fragmenting (open squares) 
and fragmenting (solid triangles) simulations carried out using 250,000 par- 
ticles and with /3sph = 2.0. The solid line, included by eye, shows a divid- 
ing line between the fragmenting and non-fragmenting cases and the grey 
region is where fragmentation can take place. At high viscosities, the dis- 
sipation is higher (Figure [TT] left panel), resulting in a faster cooling, i.e. 
a lower value of /3, required to overcome the dissipation and cause frag- 
mentation. As ogpH is decreased, the dissipation also decreases requiring 
a slower cooling for fragmentation. At very low values of ctsPH. additional 
dissipation occurs resulting in a lower value of /3 cr ; t . This may be due to ad- 
ditional particle oscillation as the shocks are not modelled adequately with 
such a low value of cksph- 



all dissipation by much (an effect also noted by |Price & Federrath| 
|20T0l >. 

In any case, Figure [5] shows that the amount of /3sph that is 
required to deal with the particle oscillations in this problem is ~ 2. 
This value ensures that the dissipation resulting from artificial vis- 
cosity is as low as possible (since the value of /3 CI it that results is 
higher) and therefore is likely to give a result that is "closer to the 
real answer". While these simulations are carried out at a single 
resolution (250,000 particles), [Bate| ( |1995) shows that ,3 S ph ~ 2 is 
sufficient to stop particle interpenetration for Mach numbers across 
a shock of ^£ « 3. In our simulations we find that the Mach num- 
bers across the shock are up to ~ 3. Note that the simulations pre- 
sented by |Rice et al.| ( |2005^ and |Meru & Bate] pOlla) were all 
carried out using /?sph = 0.2. Consequently, the calculations in 
both papers will have been affected by this and thus the converged 
value of /3 cr it is expected to be even higher than that suggested by 
Figure[T] 

5.3.2 The effect ofaspn on the critical cooling timescale 

In Section |3~T] we show that the dissipation due to the linear term in 
the artificial viscosity may play a part in the fragmentation results, 
and hence the value of /3 cr it at any one resolution. In Section |5~T| 
we show that this may indeed have been the case by considering 
the rate of convergence with increasing resolution. In this section 



we maintain a fixed resolution using 250,000 particles (i.e. keep the 
value of h/H constant) and use a fixed value of /3sph = 2.0, but 
vary the value of osph to confirm that equation|6]does indeed play 
a part in determining the fragmentation boundary. Table [3] (mid- 
dle panel) and Figure [6] summarise the simulations performed and 
the key fragmentation results. At higher values of osph the dissi- 
pation due to the artificial viscosity is expected to increase. Con- 
sequently, the cooling required to overcome this additional dissi- 
pation is larger and as a result, the critical cooling timescale for 
fragmentation is lower. As the cvsph term is decreased, the amount 
of dissipation also decreases and thus the cooling does not have to 
be so rapid, resulting in a higher critical cooling timescale. At val- 
ues below qsph ~ 0.1, however, the dissipation increases once 
again as there is not enough artificial viscosity to remove the oscil- 
lations at shock fronts. Examining the velocity dispersion of par- 
ticles in the disc around their expected almost Keplerian values, 
we find that with very low viscosity, the velocity dispersion of the 
particles increases. The particles are 'jostled' by one another when 
the viscosity is lower and the relative motions grow larger (also see 
Section [5. 3. 5} . Therefore, even though the value of q:sph is de- 
creased, the dissipation increases as the small amount of viscosity 
that is present tries to damp these larger velocities. This suggests 
that asPH ~ 0.1 is a happy medium whereby it minimises the dis- 
sipation and avoids large oscillations at shock fronts. 
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Figure 7. Graph of f3 against resolution for non-fragmenting (open squares) 
and fragmenting (solid triangles) simulations carried out with (ctsPH, 
/3 S ph) = (0.1, 2.0) using 31,250, 250,000 and 2 million particles. The solid 
line, obtained by fitting equation [T5] shows a dividing line between the frag- 
menting and non-fragmenting cases and the grey region is where fragmen- 
tation can take place. The region to the right of « 2 million particles has 
not been shaded in as it is unclear from these results alone what the shape of 
the dividing line would lie. The dotted line shows the fit using equation[l5] 
when f is set to the minimum value it can be, i.e. unity. 



Figure 8. Graph of /3 C rit against resolution of the SPH simulations carried 
out using /?sph = 0-2 (squares) and /3sph = 2.0 (triangles). The value of 
"sph is 0.1. It can be seen that the effect of increasing /3sph to 2.0 (i.e. to 
a value that minimises the additional dissipation) is to increase the critical 
cooling timescale. 



5.3.3 Determining the fragmentation boundary using optimum 
values o/aspH and /3sph 

The SPH artificial viscosity parameters that appear to produce a 
minimum excess dissipation for this problem are (qsph, /3sph) ~ 
(0.1, 2.0). However, given that the previous simulations did not use 
these optimum values (Rice et al. 2005 , Meru & Bate|201 lap , it is 
important to correct for this. We therefore carry out a number of 
SPH simulations using 31,250, 250,000 and 2 million particles to 
determine what the critical cooling timescale is for the discs simu- 
lated using (aspH, /3sph) = (0.1, 2.0). Table [3] (bottom panel) and 
Figure [7] summarise the results of these simulations. It can imme- 
diately be seen that the critical cooling timescale is higher than the 
equivalent simulations with /3sph = 0.2 (also see Figure[8j. Figure 
[9] (left panel) shows an image of a fragmented disc modelled using 
2 million SPH particles, (osph, /3sph) = (0.1, 0.2) and a cooling 
timescale of /3 = 9 which fails to fragment. Figure [9] (right panel) 
shows the equivalent disc modelled using (o?sph, /3sph) = (0.1, 2.0) 
which fragments even though it is modelled using a slower cool- 
ing time of f3 = 12. The results (Figures [7] and [8]( still show that 
as the resolution increases, the critical cooling timescale increases, 
consistent with the results presented with a lower value of /3sph- 
However, since there are only three data points with /3sph = 2.0, 
it is firstly not clear whether convergence exists and secondly, what 
function should be used to fit this data. If we assume a functional 
form as given by equation [23] we find that /3 cr it = 36.6 ± 6.9, 
assuming a first-order convergence rate (i.e. a — 1.0) as indicated 
in Section [5~T] Using the same convergence rate as found in Sec- 
tion [5TT] (i.e. a = 1.08) we find that /3 crit = 29.2 ± 2.4. This 
implies that the true value may well be as high as ~ 30. We note 



that a fit assuming a second-order convergence rate gives a poor 
fit to the data. In Section [BTT] omitting the 16 million particle data 
point makes no significant difference to the value of /3 cr i t obtained 
due to the slow convergence rate. Therefore we do not expect the 
absence of the 16 million particle data point here to significantly 
affect the fit. 

Furthermore, we attempt to fit the analytical formula from 
equation[l5] Allowing all three parameters to vary (aGi.crit, V and 
we find that aGi.crit = 0.015, Tj = 14.0 and £ = 0.3. There 
are several points to note here. Firstly, this value of the critical 
gravitational stress is equivalent to a critical cooling timescale in 
the limit of infinite resolution of /? cr it ~ 27. Secondly, the values 
of Tj and £ are much less than those obtained in Section [BTT] This 
implies that when increasing the quadratic artificial viscosity term 
to /3sph = 2.0, not only does the total dissipation decrease but 
the level of excess dissipation also decreases. However, we express 
caution here: the value of £ obtained here is lower than unity which 
is not possible. This is likely to be an artefact of using three data 
points to fit three unknowns. We therefore refit the data by setting 
£ to the minimum possible value it can be, i.e. unity. In this case, 
we find that OGi,crit = 0.024 ± 0.005 and t] = 6 ± 3 (see dotted 
line in Figure[7](. Again, we emphasise that r\ is smaller than previ- 
ously obtained in Section |5~Tj suggesting that with /3sph = 2.0, the 
excess dissipation is significantly reduced. 

We note that the analytical formula in equation [TT] is depen- 
dent on two aspects. Firstly it assumes that the ratio of the smooth- 



14 



ing length to the disc scaleheight, h/H, is given by equation 
Secondly, it assumes that the dissipation due to the artificial viscos- 
ity is indeed given by equation[8] To check the analytical arguments 
presented in Section |3~T) it is important to test these two aspects. 
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Figure 9. Surface mass density rendered image of discs modelled using 2 million particles and with Qsph = 0.1. The left image shows a disc modelled with 
/3gP H = 0.2 with a cooling timescale, (5 = 9 while the right image shows a disc modelled with /3sph = 2.0 with a cooling timescale, ft = 12. The disc 
modelled using a lower amount of artificial viscosity does not fragment even though it is modelled with a faster cooling as counterintuitively, there is excess 
dissipation with a lower value of /3sph- 
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Figure 10. Graph of the azimuthally averaged and the analytically estimated (using Equation |14) r adial profile of the ratio of the smoothing length to disc 
scaleheight, h/H , for the non-fragmenting (i.e. marginally stable, Q ?s 1) discs presented in Table|3|(bottom panel) using 31,250 (left panel), 250,000 (middle 
panel) and 2 million (right panel) particles. The analytically estimated radial profile is plotted using long dashed lines while all other lines are the simulation 
results. In the outer parts of the disc where fragmentation will occur when the cooling is close to the fragmentation boundary, the azimuthally averaged 
measured values are very close to the expected values. 



5.3.4 Testing the analytical formula for h/H 

Figure [T5] shows the analytical estimate of the ratio of the smooth- 
ing length to the disc scaleheight, h/H, against the azimuthally 
averaged radial profile of h/H for the non-fragmenting discs (i.e. 
marginally stable discs where Q ~ 1) carried out in Section [5.3.3| 
(bottom section of Table [3}. It is very clear from this graph that at 
all resolutions considered, the analytical estimate of h/H is a good 
approximation in the outer parts of the disc where the fragments 
generally form. It is important to note, that the analytical formula 
assumes the initial surface mass density profile remains constant. 
However, the discrepancy in the inner regions is due to the change 



in surface mass density profile as the disc evolves into a state of 
mechanical equilibrium on a viscous timescale. The change in sur- 
face mass density profile thus changes the value of h/H (equa- 
tion [l4j. Since the viscous timescale is shorter at small radii, the 
disc evolves more rapidly there and thus the discrepancy is larger. 
However, it is important to note that for a cooling timescale close to 
the critical one, fragmentation occurs in the outer parts of th e discs 
for surface mass density profiles shallower than E oc R~ 2 jMeru 
& Bate (20 lib I) since the resolution increases with radius (equa- 
tion [T4|. The outer parts are where the agreement is best between 
the analytically expected and azimuthally averaged values of h/H. 
We also note that the agreement between the analytical formula for 
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Table 4. Table showing the simulations carried out without self-gravity to 
determine if the dissipation due to the artificial viscosity is the same as that 
expected in a shear-dominated disc. The value of aspn is changed while 
maintaining a fixed value of /?sph = 2.0 (top panel). The bottom panel 
shows the simulations carried out with «sph = 0.1 and varying the value 
of /3sph- 

h/H and the simulation data is better with increasing resolution 
as the viscosity decreases (equations [6] and |7J and so the effective 
viscous time is larger resulting in a slower evolution of the surface 
mass density profile. We therefore conclude that the reason why the 
value of r) in Section |5.3.3| is not unity cannot therefore be put down 
to a mismatch between the analytical and actual values of h/H in 
the region where fragmentation will occur. 



5.3.5 Testing the analytical formula for /3 cr j t using 
non-self-gravitating discs 

We carry out a number of simulations of discs with 250,000 par- 
ticles with the same physical parameters as described in Section [4] 
but without self-gravity. Table [4] shows a summary of these simu- 
lations. The goal of this exercise is to see what effect the change 
in the SPH artificial viscosity parameters has on the dissipation in 
a laminar disc which should only be due to shear and whether this 
is as we would expect from the analytical formulae. To do this, 
we must start from exactly the same disc. However, the effects of 
the initial conditions must also be removed as this may affect the 
amount of dissipation. Therefore, we run a disc using a cooling 
time, j3 — 20, for 1.5 ORPs using the artificial viscosity parameters 
(«sph, /Ssph) = (0.1, 2.0) (i.e. the values that we expect would 
minimise the additional dissipation). This is equivalent to m 3 or- 
bital periods at 15 au (where this analysis is done). The cooling 
time (equation|2j is also ~ 3 orbital periods at 15 au. Since the ini- 
tial evolution time and the cooling time are approximately equal, 
the disc is then settled such that the heating matches the cooling. 
We then change the artificial viscosity parameters in the disc ac- 
cording to what is shown in Tableland run the simulations for a 
short period of time (0.1 ORPs or ~ 0.2 orbital periods at 15 au) 
and measure the total dissipation rate due to the artificial viscosity 
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Figure 12. Setting up an initial disc model with purely Keplerian shear flow, 
we measure the instantaneous viscous dissipation averaged over S3 800 par- 
ticles in a thin radial extent. We then add increasing random particle veloc- 
ities to the disc setup and measure the dissipation. In this figure, we plot the 
measured values of the dissipation divided by the analytic values expected 
for purely Keplerian flow due to the asPH (solid line) and /?sph (dashed 
line) terms separately versus the magnitude of the random velocity disper- 
sion in units of the sound speed (i.e. we plot rj (solid line) and £ (dashed 
line) as defined in equation |l 1) . The excess dissipation due to small-scale 
particle velocity dispersion can be substantially higher than that produced 
by a pure shear flow. 



in the radial range 14.9 ^ R ^ 15.1 au and compare these with 
the expected dissipation due to the artificial viscosity using equa- 
tion |C6| (using the actual values of the sound speed and smoothing 
length obtained from the simulation rather than the initial values). It 
is important to note that in order to make this comparison, we only 
calculate the dissipation over a short period of time as we are com- 
paring the instantaneous expected dissipation rate with the instan- 
taneous actual (azimuthally averaged) dissipation rate. If we allow 
the discs to run for a very long time before measuring the dissipa- 
tion rate, the discs will evolve considerably and a like-for-like com- 
parison is then not possible. This subsequent evolution takes place 
over a much smaller timescale than the orbital timescale. There- 
fore, there is no time for h/H or the velocity field of the particles 
to change. Therefore any change in the disc's dissipation must be 
due to the change in artificial viscosity parameters. 

Figure [TT| (left panel) shows a graph of how the dissipation 
rate (measured and expected) changes with qsph (using a fixed 
Psph = 2.0). The expected dissipation rate is the sum of the dis- 
sipation due to the q:sph and /jsph terms. It can immediately be 
seen that at low values of «sph, the expected dissipation due to 
the /3sph term is very important (though its contribution is often 
thought to be negligible in comparison to the asm term). At high 
values of «sph the measured dissipation matches the expected val- 
ues very well. Most strikingly, the total dissipation is higher than 
the expected dissipation from the analytical formula at low values 
of qsph - The discrepancy is about a factor of two for osph < 0.1. 

Figure[TT|(right panel) shows the measured and expected dis- 
sipation rates against /?sph (using a fixed cisph = 0.1). In this case 
the total dissipation is always approximately a factor of 2 higher 
than the expected values in these non-self-gravitating calculations. 
However, we note that there is no obvious additional dissipation at 
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Figure 11. Graph of measured dissipation rate per unit mass (short dashed line) in non-self-gravitating discs against «sph (left panel, using /?sph = 2.0) 
and /3sph (right panel, using «sph = 0.1), modelled using 250,000 particles. The expected dissipation rate per unit mass due to the cegpn (dotted line) 
and /3sph (long dashed line) terms and the combined total expected dissipation due to the artificial viscosity (solid line) are also plotted (using equation |C6fr . 
The actual measured dissipation rate is higher than the analytical estimates of the dissipation due to the shear in all cases other than when ctgpjj is high. In 
addition, the dissipation due to the /?sph term is not always negligible, as is often presumed to be the case. 



small values of /?sph compared to large values, in contrast to the 
self-gravitating calculations in Figure[5]which shows a definite dif- 
ference in results between low and high /3sph values. The /?sph 
viscosity was originally introduced into SPH to stop particle inter- 
penetration at shocks in supersonic flows. Shocks are not present in 
the non-self-gravitating calculations, so particle interpenetration is 
not an issue, but shocks play a significant role in the self-gravitating 
calculations. Thus, the apparent reduction of the dissipation in the 
self-gravitating calculations when the value of /3sph is increased 
is likely to be because particle interpenetration is stopped more ef- 
fectively with a higher value of /3sph- Regardless, both panels in 
Figure [TT] clearly show that the expected contribution to the dis- 
sipation from the quadratic artificial viscosity term can be larger 
than that from the linear term. In particular, for the simulation us- 
ing (cvsph, /3sph) = (0.1,0.2), which were the values used by 
|Rice et al7| l |2005[ > and |Meru & Bate ( 201 la I as well as many others, 
the dissipation is more than three times larger than the analytically 
expected value from the qsph viscosity alone. 

The level of dissipation expected from the SPH artificial vis- 
cosity given in Appendix [C] is lower than that measured from the 
actual simulations. What is the source of the excess dissipation that 
we find? The key is that the derivation assumes that the only con- 
tribution from the artificial viscosity to the thermal energy is due 
to shear flow in a purely Keplerian disc. Any other motions will 
add to this dissipation. In a gravitationally unstable disc, we also 
expect heating from the bulk component of the artificial viscosity 
due to shocks generated in the disc. Indeed, this is the assumed 
source of heating that is supposed to allow a gravitationally unsta- 
ble disc to achieve a quasi-steady state when an imposed cooling 
timescale, /3, is applied. However, it is exactly this maximum heat- 
ing rate that the disc can provide without fragmenting that we are 
trying to measure when we try to determine /3 c rit- The fact that 
the convergence rate of /3 cr it with increasing resolution is slow 



(first order) and that increasing /3 from 0.2 to 2.0 increases the 
critical cooling timescale significantly implies that there is a third 
source of heating. Furthermore, as demonstrated earlier in this sec- 
tion, even when we measure the dissipation in a calculation with- 
out self-gravity, we still find some excess dissipation beyond what 
Appendix [C] predicts. Figure A3 of |Lodato & R ice (2004) shows 
a calculation of the Reynolds stress in a non-self-gravitating disc 
with («sph,,3sph) = (0.1,0.2). They find that the ass param- 
eter due to this is a few xlO -3 in the range ^ R ^ 25 au. 
As a check to ensure we are consistent with previous results, we 
calculate the Reynolds stress in the same way as |Lodato & Rice] 
(2004) (though we only average over O.lORPs) and also find the 
ass value to be a few x 10~ 3 over the same radial range. Further- 
more, our results are also consistent with Forgan et al. < 20 1 lj who 
also find their ass parameter to be a few x 10 -ii (in the inner parts 
of their self-gravitating discs where the effects of self-gravity are 
not very important; see their Figure 4). 

In a non-self-gravitating calculation there are essentially only 
two possible contributions to the viscous heating. The first is from 
the Keplerian shear flow. The second is any additional particle mo- 
tions. Without self-gravity, these can only come from 'random' par- 
ticle motions. It is well known that in a typical SPH simulation the 
particles 'jostle' one another, resulting in a particle velocity dis- 
persion. This velocity dispersion results from errors in the pressure 
gradients due to the finite number of particles within a smoothing 
kernel. In a compressible SPH simulation, such motions are typi- 
cally at the level of some fraction of the sound speed. In order to 
determine the effect of these motions on the dissipation in a disc we 
perform a simple toy experiment whereby for illustrative purposes, 
we introduce different amoounts of particle velocity dispersion to 
see its effects on the dissipation in the disc. Setting up a purely Ke- 
plerian disc, we compute the instantaneous average values of the 
dissipation for particles in a small radial extent, due to the qsph 
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and /3sph terms separately. Comparing these values to those ex- 
pected from the osph and /3sph terms respectively (Appendix |C), 
i.e. 

93 2 
D a = — rjaspuCshQ (24) 

and 

Dp = ^CPspuh 2 ^, (25) 

respectively. We find that as expected, 77 = ( — 1 to a high level of 
precision (~ 1 — 3 per cent when averaging over w 800 particles). 
We then experiment with adding different levels of random veloci- 
ties in addition to the underlying Keplerian motion. The results are 
displayed in Figure[l2] where the magnitude of the particle veloc- 
ity dispersion is given as a fraction of the local sound speed in the 
disc. We see that if random motions at the level of, e.g. 30 per cent 
of the sound speed are present, the dissipation increases by factors 
of r\ = 1.7 and £ = 2.2. This provides us with an explanation 
for the excess dissipation in the non-self-gravitating calculations. 
When the level of artificial viscosity is low, the velocity dispersion 
of the particles in the disc generates a non-negligible fraction of 
the dissipation. When the viscosity is high (in particular the linear 
Qsph term), this velocity dispersion is damped, and since the con- 
tribution to the dissipation from the shear flow is larger, no signifi- 
cant dissipation beyond that expected from the shear flow is found. 
We can see from Figure [TT| that osph = 0.1 is too low to effec- 
tively damp the particle velocity dispersion, resulting in dissipation 
rates that are approximately a factor of two larger than expected. 
This does not require a high level of particle velocity dispersion - 
from Figure[l2]we see that a velocity dispersion of only w 25 per 
cent of the sound speed is enough to boost the dissipation due to 
the /3sph viscosity by a factor of two. Thus, although we show in 
Sections |5 .3. 1 1 and |5 .3 .2| that the minimum dissipation is obtained 
for (osphj/Jsph) = (0.1,2.0), this minimum dissipation is still 
larger than that expected from the analytic derivation. 

In the self-gravitating disc calculations the situation is more 
complex. Here there are gravitational forces from the fluid, shocks 
and local pressure gradients in the disc which can stir up the par- 
ticles. A low level of artificial viscosity (particularly /?sph) will 
allow particle penetration in shocks and a low value of «sph will 
be ineffective at damping post-shock oscillations and other small- 
scale particle motions. If the random motions become a substantial 
fraction of the sound speed, which they may well do since we find 
the Mach numbers across a shock to be up to « 3, the factors can 
become very large (7/ = 3 — 9 and £ = 7 — 24 for random velocities 
of 60 — 100 per cent of the sound speed). Thus, if random particle 
motions are indeed also playing a part in the self-gravitating calcu- 
lations, it should be no surprise that we infer a level of dissipation 
that is well beyond that expected from equation |C6| i.e. the counter- 
intuitive nature of /?sph that leads to Figure[5] This leaves us with a 
problem with the SPH simulations. In order to obtain a level of dis- 
sipation that is close to that predicted by a purely Keplerian flow we 
can infer from Figure[TT|that we would need to use qsph <; 1 and 
from Figure [5] that we require /3sph 2. This should cut the parti- 
cle penetration at shocks, post shock oscillations, and other particle 
velocity dispersion to low levels, thus making the analytic predic- 
tions of the dissipation accurate. On the other hand, the higher level 
of viscosity would increase the dissipation generated by the shear 
flow, thus reducing the measured value of /3 cr it at a given reso- 
lution (c.f. Figure [6](. Thus, although the simulations may be bet- 
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Table 5. Table showing the simulations carried out using FARGO and the 
key fragmenting results to test what the effect of changing the amount of 
artificial viscosity has on the critical cooling timescale. The artificial vis- 
cosity coefficient, q, is defined in equation [l6| 

ter behaved, an even higher numerical resolution would be needed 
to determine the converged value of the critical cooling timescale, 
/Scrit- We discuss other options in Section|6] 

In summary, in a purely Keplerian disc, with no random parti- 
cle velocity dispersion, the dissipation is as expected from the ana- 
lytical values in equation [C6] However, in the simulations of non- 
self-gravitating discs, particularly with low values of osph, the 
dissipation is somewhat higher than expected. We attribute this to 
random particle velocity dispersion, since there is no other source 
of heating in such discs over and above the viscous heating due 
to Keplerian shear flow. In self-gravitating discs additional particle 
dispersion will be present, which may well result in the counter- 
intuitive nature of the artificial viscosity that leads to Figure|5] i.e. 
more dissipation with lower /3sph- 

5.4 The effect of the fargo artificial viscosity on the critical 
cooling timescale 

In Section [3~2] we show that the dissipation due to artificial viscosity 
present in FARGO may play a part in the critical cooling timescale. 
Tableland Figure [T3"1surnmarise the results of the simulations car- 



© 0000 RAS, MNRAS 000, 000-000 



1 6 Farzana Meru and Matthew R. Bate 




q 



Figure 13. Graph of f3 against the FARGO artificial viscosity parameter, q, of 
the non-fragmenting (open squares) and fragmenting (solid triangles) sim- 
ulations earned out using 512 and 1536 cells in the radial and azimuthal di- 
rections, respectively. The solid line, included by eye, shows a dividing line 
between the fragmenting and non-fragmenting cases and the grey region is 
where fragmentation can take place. For low values of q, the dissipation due 
to artificial viscosity is low (and is most likely dominated by the intrinsic 
numerical diffusion) resulting in fragmentation occurring with high values 
of /3. As the artificial viscosity is increased, a faster cooling is required to 
overcome the additional dissipation, resulting in lower values of c tii ■ 



ried out to investigate this. As the artificial viscosity parameter, q, 
is increased, it becomes harder for the disc to fragment due to the 
extra heating. Consequently, the critical value of /3 required to over- 
come this and allow the disc to fragment decreases. 

At lower values of q, the effect of artificial viscosity on the 
fragmentation boundary is much less obvious. Note from Table [5] 
that when the artificial viscosity parameter is set to zero, the frag- 
mentation boundary decreases to a lower value of /3 as with the SPH 
results in Figure|6]but the effect of reducing the viscosity is much 
less pronounced in FARGO than in SPH. However, the reasoning is 
likely to be different because there is no dissipation associated with 
the artificial viscosity term since it is set to zero. We note from 
Figure [13] that /3 cr it increases rapidly as q is decreased to q ~ 0.5 
and then plateaus - this is most likely because at such low values 
of the artificial viscosity the dissipation is dominated by intrinsic 
dissipation in the code. 

For fragmentation to occur, the dissipation associated with 
high values of the artificial viscosity parameter needs to be over- 
come with a faster cooling. As with the SPH results presented in 
Sections |5 .3. 1 1 and |5 .3 ,2| artificial viscosity clearly plays a part in 
whether these discs, modelled using a grid-based code, fragment or 
not. Figure[T3]suggests that a value of q ~ 0.5 may be sufficient to 
avoid any excess dissipation. 



Simulation name 


No of 
radial cells 


No of 

azimuthal cells 




Fragmented? 


7Rfik rplk-nfi 5-hpfa3 


512 


1536 


7 


Yes 


7Sfilr rplk nfi S hpfn^ S 


S 1 9 


1 

1 J JO 


<^ 


Yes 




5 1 2 


1536 


10 


Yes 


lRf\lr rplk nfl S hftad S 


512 


1 S^n 


10 5 


No 




512 


1536 




No 


^ 1m rplk nfl S hptn 1 


1024 


3072 


14 


Yes 


1 1m rplk-nfi S-hptVl ? 


1024 


3072 


1 5 


Yes 


^ lm rplk-nO 5-hpfa1l 


1024 


3072 


16 


Yes 


1 lm rplk-nfi ^-hptn l 


1024 


3072 


1 8 


Yes 


^ lm rplk-nfi 5-hpta1 3 


1024 


3072 


20 


No 


Mm rplk-nfi 5-hpti 1 1 

1 Jlll.LClla U U. J UC Id 1 1 


2048 


6144 


20 


Yes 


1 3m_cells-q0.5-beta 1 4 


2048 


6144 


22 


Yes 


13m_cells-q0.5-betal4 


2048 


6144 


24 


Yes 


1 3m_cells-q0.5-beta 1 5 


2048 


6144 


26 


No 


50m_cells-q0.5-beta24 


4096 


12288 


24 


Yes 


50m_cells-q0.5-beta26 


4096 


12288 


26 


Yes 


50m_cells-q0.5-beta28 


4096 


12288 


28 


No 


50m_cells-q0.5-beta32 


4096 


12288 


32 


No 



Table 6. Table showing the simulations carried out using FARGO and the 
key fragmenting results. The simulations are performed using the artificial 
viscosity parameter, q = 0.5. 



100 , u 




Number of cells 



Figure 14. Graph of j3 against resolution of the non-fragmenting (open 
squares) and fragmenting (solid triangles) FARGO simulations. These sim- 
ulations are carried out with an artificial viscosity parameter, q = 0.5. The 
solid line, obtained by fitting equation|23] shows a dividing line between the 
fragmenting and non-fragmenting cases and the grey region is where frag- 
mentation can take place. The graph shows clear evidence of convergence 
of results with increased resolution. The convergence rate is second order 
with spatial resolution. The dotted line (which coincides well with the solid 
line) is obtained by fitting equation [2T| 



© 0000 RAS, MNRAS 000, 000-000 



Convergence of fragmenting self-gravitating discs 17 



100 



10 - 



10 5 



10 6 10 7 
No of cells 



10 8 



Figure 15. Graph of /3 c rit against resolution of the FARGO simulations 
carried out with an artificial viscosity parameter, q = 1.41 (squares) and 
q = 0.5 (triangles). It can be seen that the effect of reducing the artificial 
viscosity parameter to q = 0.5 (i.e. to a value that minimises the additional 
dissipation) is to increase the critical cooling timescale, with the effect be- 
ing much greater at lower resolution. 




Figure 16. Surface mass density rendered image of a disc modelled using 
FARGO with 50 million grid cells and using q = 0.5. The disc is modelled 
with a cooling timescale as high as /3 = 26 and still fragments. 



Determining the fragmentation boundary using the 
optimum value of the FARGO artificial viscosity parameter 



5.4.1 



In Section |5.2| we show that convergence appears to be reached 
at higher resolution with FARGO and that the convergence is sec- 
ond order in spatial resolution. However, these simulations use a 
value of the artificial viscosity parameter, q — 1.41, which we 
show in Section \5A\ does not minimise the dissipation. We carry 
out simulations of self-gravitating discs using a value of q = 0.5 at 



various different resolutions. Tableland Figure [14] summarise the 
simulations carried out to investigate this, and the key fragment- 
ing results. It can be seen that the effect of using a lower value 
of q is that /3 cr it is higher than obtained in Section [5^2] Figure [16| 
shows a surface mass density rendered image of one of the highest 
resolution discs (modelled using 50 million grid cells) and shows 
clear fragmentation with a cooling time as high as j3 — 26. How- 
ever, despite the critical cooling time being larger, we can see from 
Figure [14] that convergence is still being achieved. We firstly fit 
the data using equation [23] and find that /3 cr it = 28.0 ± 0.2 and 
a = 1.89 ± 0.05 showing that the convergence rate is second order 
with spatial resolution. We then fit the data using equation [27] and 
find that aoi.crit = 0.0145 ± 0.0001 and £ = 3.16 ± 0.04. In the 
limit of infinite resolution, this value of aGi,crit is equivalent to a 
critical cooling timescale, /3 cr it « 28 (using equation[4](. 

Figure[l5]shows the fragmentation boundary (with error bars) 
using q = 1.41 (as in Section [5T2] > and q = 0.5. As the resolution 
increases, the difference between the two sets of results decreases: 
since the convergence with FARGO is fast, i.e. second-order, the ef- 
fect of using different values of the artificial viscosity parameter 
(i.e. q = 0.5 versus q — 1.41) becomes negligible with 50 million 
grid cells compared to a lower resolution. This further corroborates 
that at a higher resolution, the artificial viscosity plays less of a 
part in the fragmentation boundary. As suggested by the analytics 
in Section [3~2] it is expected that the artificial viscous dissipation 
should decrease both when the resolution is increased and when 
q is decreased (until the numerical dissipation becomes dominated 
by intrinsic grid dissipation). Indeed, the higher value of £ obtained 
here in comparison to that in Section |5.2| suggests that the dissi- 
pation due to the artificial viscosity has been minimised and that 
the intrinsic grid dissipation is becoming more important, consis- 
tent with Figure [T3] We note that the value q — 0.5 effectively 
means that the shock is spread over approximately half a grid cell 
which will affect the treatment of shocks. We emphasise that we 
choose this value since it gives the lowest artificial heating rate that 
is possible with FARGO, as done so with the SPH simulations. More 
importantly, we show that the choice of the value of q has much less 
of an effect at higher resolution than at low resolution. 



6 DISCUSSION 

The non-convergence of results concerning the fragmentation of 
self-gravitating discs has opened up a number of questions con- 
cerning both the physics and the numerics involved in determin- 
ing whether a disc will fragment into bound objects. Consequently, 
|Meru & Bate ( 201 la) presented the possibility that either the criti- 
cal cooling timescale was larger than originally thought, or the ex- 
treme possibility that the physics behind the fragmentation of discs 
needs to be reconsidered. In this paper we find that using both SPH 
and the FARGO codes i.e. a three-dimensional particle-based La- 
grangian code and a two-dimensional grid-based Eulerian code, re- 
spectively, the artificial viscosity that is used to accurately model 
shocks plays a significant part in the convergence rate. Not only 
do we show that convergence can occur, but we also show that the 
rate at which it occurs is as expected from analytical arguments 
involving artificial viscosity (first-order with linear resolution for 
SPH and second order for FARGO i.e. a faster convergence with 
FARGO). This affects the results on the fragmentation boundary in 
both SPH and grid-based calculations. 

In particular, we conclude that oscillations at the shock front 
and particle interpenetration may not have been adequately ac- 
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counted for in the previous SPH simulations (via the quadratic ar- 
tificial viscosity term, /3sph). Those SPH simulations that used a 
value of /Ssph that was too low, may have counterintuitively re- 
sulted in more dissipation, causing fragmentation to have been un- 
derestimated. After minimising the additional dissipation associ- 
ated with the artificial viscosity employed in both codes, we find 
that the critical cooling timescale is at least as high as 20 and per- 
haps even as much as w 30, for a ratio of specific heats, 7 = 5/3. 

Previous simulations that investigate the effects of artificial 
viscosity show that the strength of the gravitational instabilities 
are weakened ( [Pickett et al. 2000 1 and clump formation is reduced 
l |Boss|20 06:> when artificial viscosity is employed. These results are 
in the same sense as our results, i.e. excess effective viscosity re- 
duces the propensity for fragmentation. 

|Mayer et al.| f2004 1 carried out a test on the fragmentation and 
disc evolution resulting from the inclusion of artificial viscosity in 
their three-dimensional SPH code. Pickett & Durisen (2007) car- 
ried out a similar test using a three-dimensional grid-based code. 
Both sets of authors perform their tests on isothermal simulations 
i.e. they only include the effects of artificial viscosity in the mo- 
mentum equation and did not consider its heating effects in the en- 
ergy equation. They suggest that artificial viscosity may reduce or 
even prevent clump formation from occurring. While our results are 
consistent with theirs with respect to preventing clump formation 
when artificial viscosity is increased, we stress that the dissipation 
associated with the artificial viscosity plays a key role in the frag- 
mentation results. 

Another possible numerical parameter that may affect the re- 
sults is the gravitational softening used in the two-dimensional 
grid code. Mull er et al.| ( |2012| l show that an incorrect value of 
the gravitational softening length in two-dimensional disc simu- 
lations can significantly affect the fragmentation conclusions: a 
low value causes the gravitational forces on short distances to be 
over-estimated, resulting in the conclusion that fragmentation does 
occur, when the converse conclusion is reached for larger values 
of the softening parameter. Indeed, it is well known that three- 
dimensional discs are more stable than two-dimensional discs since 
the vertical component dilutes the effect of gravity (Toomre 19641. 
Thus, incorrectly taking into account the effects of the vertical di- 
rection in a two-dimensional simulation may cause the disc to be 
more prone to fragmentation than its three-dimensional equivalent. 
Muller et al. (2012) show that a value of « Q.QH is required to 
model the gravitational forces correctly (though they do say that 
a comparison with 3D simulations is required). Since we are us- 
ing a softening length of 3 x 10~ 4 H, our FARGO simulations may 
overestimate fragmentation. 

For the SPH simulations, although we see evidence for con- 
vergence of the critical cooling timescale, the convergence rate is 
only first order with increasing resolution. This is partly due to the 
larger dissipation than that predicted by the continuum limit of the 
SPH equations in a shear flow. We argue that the excess dissipa- 
tion is due to small-scale particle velocity dispersion. In non-self- 
gravitating discs, this results from pressure gradient errors due to 
discretisation, but in the self-gravitating discs there are other po- 
tential sources: primarily particle penetration at shock fronts and 
post-shock oscillations (particularly when the levels of artificial 
viscosity are low), but perhaps also discretisation errors in the self- 
gravity. To achieve well-behaved dissipation (i.e. that which is close 
to that predicted by the continuum limit of the SPH equations) re- 
quires asPH 1 and /3sph <; 2. However, the dissipation from the 
shear flow is then relatively high meaning that even higher resolu- 



tion would be necessary to obtain a converged value of the critical 
cooling timescale, /3 cr it. 

However, there are many possibilities that might improve the 
SPH performance. We have employed the most basic form of SPH 
artificial viscosity in this paper (i.e. constant values of osph and 
/Ssph). An obvious aspect to investigate is whether a viscosity 
switch such as those proposed by Morris & Monaghan ( 1997) or 
ICullen & Dehnenl polo) which increase qsph and /3sph in the 
presence of a shock and allow them to decay away from shocks can 
improve the convergence rate. This could potentially provide high 
viscosity to avoid particle penetration at shocks and post-shock os- 
cillations (i.e. reducing small-scale particle velocity dispersion), 
but retain low viscosity (and thus low heating rates) in the bulk 
of the disc, thus minimising the heating due to the shear flow. Since 
some of the particle velocity dispersion originates from pressure 
gradient errors, another possibility is to try a more accurate kernel. 
The quintic spline kernel generally performs better than the cubic 
spline kernel and Morris et aL] J1997) reported that it significantly 
reduced velocity field noise in their calculations. Testing these vari- 
ants is beyond the scope of this paper, but we expect that these and 
other SPH variations may be able to significantly improve the per- 
formance of SPH on this problem. We stress that these possible 
improvements to the SPH convergence rate will not decrease the 
value of the critical cooling timescale and thus the values obtained 
in this paper indicate lower limits for /3 C iit- 

Since [Meru & Bate (2011a) published their results highlight- 
ing the convergence problem, a number of authors have attempted 
to explain the non-convergence. Lodato & Clarke (2011) specu- 
lated that the cause may be the artificial smoothing of the density 
enhancements in SPH or a larger than expected level of artificial 
viscosity. Our results clearly show that artificial viscosity plays a 
major role in numerical determinations of the critical cooling rate. 

|Paardekooper et aT] ( |2011[ l suggested that the boundary be- 
tween the turbulent inner disc region and the laminar outer disc 
region (a natural consequence of starting with smooth initial con- 
ditions) may cause an edge in the disc that becomes more and 
more pronounced at higher resolutions, making it easier to frag- 
ment. They suggested that if the smooth initial conditions were re- 
moved, convergence could be achieved. In light of the new results 
presented in this paper, the effect of edges should be considered 
in more detail. If edge effects do play a part, it is unclear whether 
they should always continue to become sharper at higher resolution 
(and hence inconsistent with the results presented here), or whether 
they should "saturate" at higher resolution (and thus consistent with 
these results). It is important to note, however, that Bate (201 1 1 per- 
formed radiative transfer calculations of molecular cloud collapses 
and found that disc fragmentation is more prevalent in higher res- 
olution calculations. These discs did not begin with smooth initial 
conditions, and yet a similar resolution dependence was seen. 

More recently, Paardekooper (2012) carried out shearing sheet 
simulations, similar to those performed by Gammie (2001 1 where 
no such edge effects should play a part. He found that as the reso- 
lution was increased, the critical cooling timescale also increased, 
showing that the convergence issue is not restricted to global simu- 
lations, but also affects local simulations. He found fragmentation 
for at least as large as /3 cr jt ~ 20. However, this was for simu- 
lations carried out with a ratio of specific heats, 7 = 2. This is 
equivalent to a maximum gravitational stress as least as small as 
c*Gi,crit « 0.011, consistent with the value of the gravitational 
stress found using our global simulations. 

It is important to note that many of the previous simulations 
that have attempted to explain the convergence problem highlighted 
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by |Meru &"B ate (201 la) have tried to do so by carrying out simula- 
tions with resolutions in the non-convergent region of the resolution 
space shown here in this paper. It is therefore hard to interpret those 
results since they may have been affected by artificial viscosity. For 
those codes that do not use artificial viscosity, other sources of nu- 
merical diffusion relating to the specific implementation may be 
important and the effect of these need to be thoroughly understood. 
It would be interesting to try to understand the convergence prob- 
lem with a Godunov scheme that does not implement an artificial 
viscosity, or to apply a fixed Navier-Stokes viscosity. 

While this work focusses on the effects on fragmentation of 
self-gravitating discs, the key conclusion that artificial viscosity can 
play a significant role in the physical interpretation of simulations is 
more general. We emphasise that any simulations whose outcome 
is highly dependent on the thermodynamics of a problem should 
ensure that the effects of artificial viscosity in their code are well 
understood as well as highlighting the importance of convergence 
of results with both resolution and numerical method. 

6.1 Implications for the fragmentation of real discs 

A critical cooling timescale of /3 cr it ~ 20 or ~ 30 is equivalent 
to a maximum gravitational stress of aGl.crit ~ 0.02 or ~ 0.013, 
respectively. Clarke (2009) produced an analytical model for the 
structure of a gravitationally unstable disc which is subject to real- 
istic cooling. She showed that for optically thick discs that are suf- 
ficiently low in temperature that they are dominated by ice grains, 

«GI = 0.4 (zr^— Y, (26) 
\ 100 au J 

for a disc with interstellar opacities and surrounding a IMq star, 
where R is the radius being considered. This relationship shows 
that for a maximum value of the gravitational stress, a critical ra- 
dius, i? cr it, can be found outside of which fragmentation can occur 
(for a disc with a shallow surface mass density profile). While the 
previously accepted result of QGi.max ~ 0.06 gives a critical ra- 
dius of i?crit ~ 68 au, the values of /3 cr it obtained here moves the 
critical radius inwards to i? cr it ~ 47 — 51 au (for a disc around a 
IMq star using interstellar opacities). The core accretion scenario 
is thought to occur out to ~ 10 au, while gravitational instability is 
historically thought to operate outside of ~ 70 — 120 au (Rafikov 
|2009| [Clarke 2009). Therefore, an intermediate radial region exists 
where no one in situ formation method adequately seems to de- 
scribe the formation of planets. Our results show that this gap can 
at least partly be bridged if the true critical cooling timescale is as 
much as /3 C rlt ~ 20 — 30. 

We point out that equation [4] is derived by assuming that the 
dominant form of heating in a disc is that due to the gravitational 
instabilities. In a real disc, there may be a contribution to the stress 
from the magnetorotational instability (MRI), omri- In this case, 
equation|4]may be written as 



9 7(7 — 1) (aGI + "MRl) ' 

Therefore, while the critical cooling timescale for a purely gravita- 
tionally unstable disc is quite large, if the contribution to the grav- 
itational stress from the MRI (or in fact other heating sources) be- 
comes important, a faster cooling will be required to overcome this 
additional heating and allow the disc to fragment. Since we find that 
the critical stress may be as low as aGi,crit ~ 0.01, the heating due 



to MRI will certainly be expected to be important if it provides an 
effective stress of approximately this level or higher. Even if ckmri 
is a factor of 10 smaller, this will still make a 10% difference to the 
heating which can change the critical cooling timescale required 
for fragmentation. 



7 CONCLUSIONS 

We perform hydrodynamical simulations using a three-dimensional 
Smoothed Particle Hydrodynamics code and a two-dimensional 
Eulerian grid-based code of self-gravitating discs to investigate 
how the presence of artificial viscosity may affect fragmentation 
results. We present additional SPH results to those presented by 
|Meru & Bate (2011a) as well as perform similar simulations us- 
ing the grid-based hydrodynamics code, FARGO, and show that 
convergence with resolution of the critical cooling timescale can 
be achieved with both codes. We show that the previous non- 
convergent results are largely due to the effects of artificial viscos- 
ity that play a more prominent role at lower resolution. We find that 
the convergence rate of the critical cooling timescale required for 
fragmentation is first order in spatial resolution using SPH and sec- 
ond order using FARGO. Furthermore, we find that the dissipation 
from the artificial viscosities in SPH is exactly as we would ex- 
pect in a purely laminar disc. However, if random particle motions 
are present, they can produce dissipation due to artificial viscosity 
that is larger than expected. In self-gravitating discs, using a value 
of the quadratic artificial viscosity term that is too low can result 
in counterintuitively high dissipation. This may be caused by ad- 
ditional random particle velocity dispersion due to the presence of 
shocks, causing the dissipation to significantly deviate from that 
expected from the SPH continuum limit equations. In addition, the 
dissipation due to the /?sph term may not be as small as previously 
assumed and should not be ignored. 

We show using analytical arguments and numerical simula- 
tions that as the resolution is increased, the artificial viscosity term 
becomes less important and the rate of convergence is as expected 
from the analytical arguments. With the particular setup adopted 
here, which is the same as that used by Rice et al. (2005} and |Meru] 
|& Bate| ( |2011a| ), we find that once the effects of artificial viscosity 
have been minimised, the critical cooling timescale converges with 
increasing resolution to a value at least as high as /3 cr it ~ 20 and 
perhaps even as high as /3 cr it ~ 30. However, a convergence be- 
tween the two codes has not yet been achieved. We conclude that 
this is much more of a problem than had previously been supposed, 
in part due to the slow convergence rate of SPH and in part due to 
the enormous resolution required to obtain convergence. The criti- 
cal cooling timescale is a factor of « 3 — 5 times larger than the 
value of /3crit ~ 6 that has been used in the past, and is equivalent 
to a maximum gravitational stress of QfGi.crit ~ 0.013 — 0.02 that 
a disc can handle before it fragments (in contrast to the previously 
obtained value of QfGi.crit ~ 0.06). 

We show that using values of the artificial viscosity that do not 
minimise the additional dissipation caused by it can significantly 
affect fragmentation results and we expect that any other results 
that sensitively depend on the thermodynamics of a problem, e.g. 
collapse of AGN discs and molecular clouds into stars may also be 
affected. This highlights the importance of ensuring that the artifi- 
cial viscosity does not play a significant role when carrying out nu- 
merical simulations. We show that fragmentation of self-gravitating 
discs can be suppressed if the effects of artificial viscosity are not 
carefully considered. This suggests that fragmentation of discs into 
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bound objects (e.g. for the formation of planets, binary companions 
and stars formation in galaxy simulations) is easier than previously 
thought. 
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APPENDIX A: ANALYTIC DERIVATION OF THE SHEAR VISCOSITY PRESENT IN SPH 

In the following sections, for the sake of clarity, we make the simplifying assumptions that the SPH particle smoothing length h, the sound 
speed, c s , and the density, p, are all slowly varying (i.e. are constant). 



Al The original SPH artificial viscosity 

The standard SPH artificial viscosity method described by Monaghan & Gingold ( 1983 ), which is a time-independent fixed artificial viscosity, 
adds the following term to the momentum equation: 

' iVi = -^mjHijViWij (Al) 



dt 



where 




(A2) 



m,j is the mass of particle j, Wij is the smoothing kernel adopted, h is the smoothing length and vy = Vj — Vj is the velocity difference 
between particles i and j. The quantity rj 2 — O.Ol/i 2 is included to avoid divergence for small separations between neighbouring particles. 
The viscosity involves two terms, the strengths of which are controlled by the parameters qsph and /?sph- Note that the artificial viscosity 
is only applied when particles approach each other and is turned off when they recede from each other. 

It has been recognised for some time that in the continuum limit, the qsph viscosity term applies both a bulk and a shear viscosity 
which has the form of a Navier-Stokes type viscosity (Monaghan 1985, Pongracic 1988; Meglicki et al. 1993). In particular, Meglicki et al. 
{1993} provide a clear derivation showing that the viscous acceleration due to the qsph viscosity is given by 

= *p [ ' {CsP } + V ( Cs <° V ' v )l ' (A4) 

where 

b%1 ~dxi + dx^ (A5) 

is the deformation tensor. The first term in equation |A4| is a shear viscosity, while the second term is a bulk viscosity. The constant, re, depends 
on the number of spatial dimensions and the kernel used by the SPH code. If the qsph viscosity is applied in three dimensions and between 
both approaching and receding particles (unlike equation|A2[l then (see Appendix|Bl[> 



4n f sdVK .... 

-is r^ dr - (A6) 



The value of this integral depends on the kernel that is being used. We use the standard cubic spline kernel 

1 - fg 2 + fg 3 forO sC g < 1 

W(q,h)= ° 



h d 



|(2 -g) 3 forlsCg<2, (A7) 

otherwise, 



where d is the number of dimensions, a is the normalisation constant equal to 2/3, 10/ (7n) and 1/n in one, two and three dimensions 
respectively, and q — r/h. In this case, the integral has the value — 3/(47r), such that n — 1/5. We find that in general the shear viscosity 
contribution to the momentum equation in two and three dimensions can be written 

^=,V.S=^qsphc s W.S, (A8) 

where v is the kinematic shear viscosity. Thus, for exampl e, when simulating acc retion discs, |Artymowicz & Lubow| ( [1994) used 



|qsphc s /i for their two dimensional SPH simulations, while |Lodato & Price 1 2010 1 used v = ^qsphc s /i for their three dimensional SPH 
simulations. 

However, as expressed in equation |A2| it is usual in SPH simulations to only apply the artificial viscous force between approaching 
particles. Thus, in general the kinematic shear viscosity in SPH simulations is a factor of two smaller and, for three dimensional calculations, 
is 

v = t^qsphc s /i. (A9) 
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A2 A more recent variation of SPH artificial viscosity 

Recently, a slightly different form of artificial viscosity has been applied in SPH codes (Chow & Monaghan 1997 , Price & Monaghan 2004) 
where, essentially, equation |A3| is replaced by 

(Hj = Vij ■ iij. (A10) 

This is the form of the artificial viscosity that we use and is constructed in analogy to dissipative terms in Riemann methods. This form also 
avoids the arbitrary quantity 77 2 which was included to avoid numerical divergences. It also means that the magnitude of puj differs by a 
factor of Vij /h from the original form of the viscosity. Note that different implementations can differ slightly in the way that the coefficients 
qsph and /Ssph enter the equations. For example, in some implementations the coefficient of the quadratic viscosity term is actually given 
by the product of a and /3, while in others the coefficient is fixed to be 2a. Some implementations may also differ in the value of Hij by 
factors of two, so care needs to be taken when evaluating the continuum limit of the viscosity in any particular SPH code. 

It can be shown (see Appendix |B}, that when using equations | A2| and [A 1 0| and the standard cubic spline kernel in three dimensions that 
the shear viscosity is actually 18% larger than for the original artificial viscosity such that 

31 

v = ^^aspHCsh. (All) 



A3 The /3sph viscosity 

Although many past studies have considered the continuum limit of the linear osph viscosity in SPH, to our knowledge, nobody has 
considered the contribution of the quadratic /?sph viscosity. By inspection of equations |A2| and [A3] we can determine how the shear viscosities 
due to both the qsph and /Ssph should scale. For a pure shear flow, p,ij provides an estimate of the shear rate of the fluid multiplied by h. 
Thus, since the kinematic viscosity is the ratio of the shear stress (given by pHij) to the shear rate, we expect the kinematic viscosity due to 
the qsph term to scale as 

u a = q a a S pHC s h (A12) 

where q a is a constant of proportionality. This is consistent with the analysis given above. An alternative method of arriving at this equation 
is to use the fact that from kinetic theory, kinematic viscosity is proportional to the characteristic speed of interchange between particles and 
to the characteristic distance over which interchange occurs. In this case, the characteristic speed of particle interchange is given by c s and 
the distance over which particles interact is a smoothing length, h. 

The | von Neumann & Richtmyer ( 1950) type viscosity (SPH /3-viscosity) is a second order viscosity with the viscous forces depending 
on the square of the relative speed of the particles. Therefore, although the characteristic distance over which the viscosity acts is still a 
smoothing length, the characteristic speed is now the relative speed between particles over a smoothing length. In a flow directed in the 
x-direction which is sheared in the y-direction this is given by 




(A 13) 



and thus the /3-viscosity is expected to scale as 



vp = qpPsmh 2 [ ^ ) (A14) 



dv x 



where qp is a constant of proportionality. 



APPENDIX B: SPH ARTIFICIAL VISCOSITY IN THE CONTINUUM LIMIT 

One way to evaluate the constants of proportionality, q a and qp, is to use the defining equation for kinematic viscosity. This can be obtained 
by considering the shearing force produced by a viscous fluid on a plane running parallel to the direction of motion of the fluid. If the fluid 
flows in the x-direction and there is a velocity gradient across the flow in the y-direction, then the kinematic viscosity of the fluid is defined 
by 



A= vp ^ (B1) 

which gives the force per unit area exerted on the plane surface as a function of the kinematic viscosity, v, the density of the fluid, p, and the 
shear rate of the fluid. 

The force exerted on a volume element of fluid is determined by considering two planes parallel to the flow. If the two planes are parallel 
to the x-z plane and they are separated by a distance 5y, then the net force on the fluid element is 
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Fx - Fa = 



fipi 



SxSySz 



(B2) 



so that the force per unit volume is given by 



F 

r 



d 2 v x 
dy 2 



(B3) 



in the limit that Sz — > 0, assuming the density and kinematic viscosity are constant, which is a valid assumption to make if the two planes 
are sufficiently close to each other. The force per unit volume for a fluid element is simply the acceleration of the fluid element multiplied by 
its density. Thus, the equation 



d v- 



dt V dy 2 



(B4) 



is obtained, which can be compared directly to the SPH momentum equation |Al| Using equations | A 1 2| and | A 1 4| we therefore find that the 
specific force can be expressed as 



dvx 
dt 



qc,aspuCah 



dV 
dy 2 



and 



(B5) 



dvx 



/3SPH ^ " 



d 2 v x 
dy 2 



for the qsph and /?sph terms, respectively. 



(B6) 



Bl Evaluating the constant of proportionality, q a , with the original form of the SPH artificial viscosity 



To evaluate the constant of proportionality for the aspH term in the artificial viscosity, we can derive the continuum limit and compare it to 
equation|B5|to determine the magnitude of q a . Using equations|Al|-[A3| the force per unit mass due to artificial viscosity on particle i is 



dvi 

~d7 



Vjj ■ Tjj 2 (Vij ■ Vij) 

~ctspunc s — ■ — — + Psph/i 



QSPH^C S 



2 {Vij -Vijf 



' Kernel 

The artificial viscosity in the continuum limit due to the linear term (qsph) is. 

dvi 

di " ./Kernel (if, ' + V 2 ) 



(4) 2 



Vl3 ' r ' J VMnd : \c. 



2 \ v * " »J l 



Following Appendix A of Meglicki et al. ( 1993), we expand Vj around to give 



p dvi &x v Ax q d 2 Vi Ax p Ax q Ax a d 3 Vi 
X dxP + 2 dxPdxi + 6 dxPdxidx a 



In addition, 



dW zj Ax k 
dr r 



(B7) 



(B8) 



(B9) 



(BIO) 



where r = \rij\ = \ri — rj\. Inserting equation BIO into the first part of equation B8 and retaining terms of order 0(h ), which are the 
lowest non-vanishing terms, gives 



dvi 
dt 



aspHhc s 



dv5 1 . „ 9 »[ 

L j Aa; - 

8xp 2 dxPdxi 



A^A^A^^d 3 ^ + 0(h 6 ). 
dr 



(Bll) 



We note that integrating a term with (Ax)', where t is odd, over a symmetric kernel yields a zero result. Therefore, the result of the integration 
of the first term in equation |Bl l| is zero. Simplifying equation |B 1 1 | gives 



dvi 
dt 



aspuhcs d 2 v r j 
2 dxPdxi 



Ax p Aj; r Ax fc A-' J 



1 dWy j3 



r 3 dr 



dx. 



(B12) 
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The integral is a fourth order symmetric isotropic tensor which can be written in the form 



-f 

J Hi 



Ax p Ax r Ax Ax' 



1 dWi 



ij A 3 



dr 



d X = K a (S pq S r k + Sp r S q k + SpkSrq) 



where 



A 



dr 



dr. 



Equation|B12|can therefore be written 



dv[ 
dt 



KaCtSPHtlCs 



hc B d 2 Vj 



dxPdxi 



(fipq5 r k + SprSqk + SpkSrq)- 



Contracting with the delta terms yields 



dvj ^ K a aspnhc s 
dt ~ 2 



2„P 



0V a^: 
■ + 



+ 



d 2 v\ 



dxPdxP dxPdx k dx k dx r 



(B13) 



(B14) 



(B15) 



(B16) 



If we now assume that locally there is a constant flow in one direction (e.g. the x-direction) with a velocity gradient in an orthogonal direction 
(e.g. the y-direction) so that dvf / dx k = and d/dx k (dv k / dx p ) = 0, as we would expect for a shear flow, then the above equation can be 
simplified to give 



dvl 
dt 



t a Q S PHCsft / d 2 Vj 

2 V dxPdxP 



(B17) 



Comparing this to equation B5 yields K a /2 = q a 



Equation |B 1 3 1 defines n a . To calculate K a in three dimensions, we need to sum over all possible combinations of r, p and q in equa- 
tion |B13| The simplest case that yields a non-zero value of the right hand side of equation [B 1 3 | involves k — r = p — q such that equation [B13| 

gives 



(Ax K ) 



dr 



(B18) 



' Kernel 

There are three additional cases where k is equivalent to one other letter while the remaining two letters are equal but in an orthogonal 



direction to k. Since in three dimensions there are two orthogonal directions to k, equation|B13|gives 



(Ax k ) 2 (Ax 



q\2 



1 dW, 
r 3 dr 



H ,3 



d x = 6n a . 



Summing equations |B18| and |B19| together yields 



' Kernel 

Without loss of generality, we use x = r sin 6 cos 



[(Ax fc ) 4 + 6(A 3 ; fc ) 2 (Aa ; 9 ) 2 )] 



1 <UVi <uV ■),,.. 



r 3 dr 



Ax , y = rsin^sini 



Ax 9 , and d 3 x = r 2 sin6» dr dO 



(B19) 



(B20) 



and substitute 



into equation B20 For the ^-component we integrate over 6 = [0, n]. However, for the (^-component, care must be taken to integrate over 



the correct range since we only consider particles that are approaching each other. Figure |CT| shows that in the frame of the particle being 
considered the ranges <j> — [0, 7r/2] and (f> — [it, (3n)/2] involve particles approaching each other whereas outside this range particles recede 
from each other. Integrating yields 



2tt 
15 



' Kernel 

Using the standard cubic spline kernel given above (for three dimensions) 

_ a dW(q,h) 



dr 



3 dW 

r dr. 

dr 



3^ 

47T ' 



Therefore, substituting into equation |B2 1 | yields K a = 1/10, and so equation |B 16| becomes 

dv x 
dt 



1 , d 2 v x 

_ aspHC ^_ 



(B21) 



(B22) 



(B23) 



We note that the same constant is achieved if the integration in the (^-direction is done over all space and simply divided by two. 



B2 Evaluating the constant of proportionality, q a , with the recent form of the SPH artificial viscosity 



For the more recent form of SPH artificial viscosity (equation |A10[ >, the procedure is identical, but the expansion and, thus, the integral is 
slightly different. The required expansion is 

a 3 ^ 



pdv^ Ax p Ax" d 2 Vj Ax p Ax q Ax' 
X dxP + 2 dxPdxi + 6 dxPdxidx a 



Ax r 



(B24) 
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central star 



= 7r/2 




= 



4> = 3tt/2 

y 



>• X 

Figure CI. Schematic diagram showing the fluid motion (denoted by solid arrows) surrounding an SPH particle (black dot) in its rest frame, as the gas disc 
orbits the central star. In the region of </>-space where <j> = [0, 7r/2] and <j> = [n, (3tt)/2] (shaded region) the fluid is approaching the SPH particle and in 
this region artificial viscosity is turned on while in all other areas the fluid is receding from the SPH particle and here the artificial viscosity is turned off (see 
equation |A2( . The origin is the location of the SPH particle. 



so that equation |B 12| becomes 

dwf - QSPHCs Q2< f A/A/A/A.' ', <m •<! ,B25> 



dt 2 dxPdxl i K ernel ' ?* 2 dr 

Note that h is missing from this equation and the integral differs by a factor of r. Thus, the equivalent of equation |B21| that needs to be solved 
is 

idWi, 



'I 



crncl dr 



dr (B26) 



where A = (— 2tt)/15 (ensuring that we account for the fact that the viscosity is only applied between approaching particles), but this time 
the integral has a value of — 31/i/(357r). Therefore, we obtain 



dv x 
dt 



31 , d 2 v 



asPHC B h^-^-, (B27) 



a SPH 

which as noted above is approximately 18% larger than the value obtained for the original form of the artificial viscosity. 



APPENDIX C: THE DISSIPATION ASSOCIATED WITH a SPH AND /3 SPH 

For this present paper we are more interested in the magnitude of the thermal dissipation provided by the viscosity than the angular momentum 
transport as such. The dissipation rate per unit mass in a viscous accretion disc is given by 

? = ^ - - (CI) 

dt E V dRJ dR V dRJ 4 ' 

where Trj, is the stress tensor and the final equality assumes a Keplerian disc. In this case the shear rate is 

Rather than derive the dissipation rate from the continuum limit of the SPH momentum equation, we can derive the dissipation rate 
directly from the continuum limit of the SPH energy equation 

~dT = \ " • v ■ ViWij. (C3) 

i 

Taking the continuum limit of this equation and using the more recent form of the artificial viscosity gives 

-jT ~ - / [-«sphc s v • r + ^sph(v ■ r) J v • r — — d x. (C4) 

d * ^ ^Kernel dr 
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The simplest way to obtain the dissipation rate due to the artificial viscosity in a Keplerian disc is to consider a small patch of the disc 
around a particle such the the local velocity field is given by 

v x « — , (C5) 

where fio is the angular velocity of the fluid at the radius being considered and y is the displacement in the inward radial direction (see 
Figure |CT) . Inserting this into equation |C4] and taking x = r sin 9 cos <j> and y = r sin 6 sin <f>, we obtain 

du ^ 9 
dt ~ 4 



525 707T 



(C6) 



taking care to integrate only over the ranges cf> — [0, 7r/2] and = [ir, 3n/2] (i.e. where the flow is approaching and not receding; see 
Figure [CT| . 

For the original form of the artificial viscosity, the dissipation is 

du ~ 9 
dt ~ I 



20 357T 



(C7) 



again only taking the integral over the regions where the flow is approaching. 

Note that the coefficients preceding q:sph in equations |C6| and |C7| are the same as the coefficients appearing in the the kinematic 
viscosity as given by equations |A9| and [Al 1| respectively, as is expected from equation |Cl| 

Using H — Cs/Q,, the ratio of the dissipation rates associated with the linear and quadratic artificial viscosity terms using the recent 
form of the SPH artificial viscosity is given by 

D a 627r «sph H 



D 135 /3sph h 



(C8) 



Typical values of «sph and /?sph are frequently within a factor of 2 of each other (with /Ssph = 2aspH). Therefore, the important variable 
that determines the relative magnitude of the dissipation associated with the linear and quadratic SPH artificial viscosity terms is the ratio 
of the smoothing length to disc scaleheight, h/H. In a well resolved disc, h/H <C 1 such that D a 3> Dp. However, if the disc is poorly 
resolved and/or /3sph 3> «sph the dissipation associated with the quadratic viscosity may be significant. Finally, note that the above 
equations assume that the only viscous dissipation comes from the shear in a purely Keplerian disc. 
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